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ABSTRACT 

Four  techniques  for  the  numerical  solution  of  partial 
differential  equations  and  eigenvalue  problems  were 
investigated.   Typical  problems  considered  were  elliptic 
partial  differential  equations  of  the  form 

U   +  U   =  f(x,y),  (1) 

xx    yy      9J    * 

or 

U   +  U   +  X2U  =  0,  (2) 

xx    yy  » 

where  appropriate  boundary  conditions  are  specified  so 
that  the  problem  is  self-adjoint. 

The  four  methods  are  relaxation,  Galerkin,  Rayleigh- 
Ritz,  and  dynamic  programming  combined  with  Stodola's 
method,  for  eigenvalue  problems. 

The  results  indicated  that  for  eigenvalue  problems 
relaxation  or  dynamic  programming  modified  is  to  be 
preferred  usually  and  for  partial  differential  equations 
Galerkin  or  dynamic  programming  is  preferred. 
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I.   INTRODUCTION 

Initially  for  this  thesis  it  was  planned  to  investigate 
methods  for  finding  eigenfunctions  and  eigenvalues,  with 
a  particular  interest  in  the  oscillation  in  basins  such  as 
harbors  and  bays.   The  report  by  Angel  [10]  introduced 
a  new  mothod,  that  of  dynamic  programming,  for  the  solu- 
tion of  some  partial  differential  equations.   This  method 
seemed  promising  and  it  was  then  extended  here  to  the  solu- 
tion of  other  partial  differential  equations.   It  also 
made  it  possible  to  invert  a  linear  differential  operator 
and  hence  apply  Stodola's  method  in  finding  eigenvalues 
and  eigenfunctions.   This  led  naturally  to  a  comparison 
with  several  procedures  already  known  to  try  to  compare 
convergence,  speed,  and  accuracy,  by  applying  them  to  the 
solution  of  several  rather  simple  problems. 

It  is  assumed  that  the  reader  has  some  familiarity  with 
the  techniques  of  replacing  differential  equations  by  dif- 
ference equations  and  the  standard  technique  of  separating 
variables  in  linear  partial  differential  equations.   The 
regions  and  their  boundaries  were  assumed  to  be  "nice," 
not  excessively  irregular.   Some  knowledge  of  calculus 
of  variations  also  is  desirable  but  not  necessary;   ref- 
erence (5}  provides  more  than  adequate  background. 

Four  problems  were  divided  into  two  categories.   In  the 
first  category  were  eigenvalue  problems.   Two  typical  ones 
were 


u   +  u   =  -  X2U 
xx    yy 


subject  to  the  constraint 


in  A, 


(1.1) 


U(x,y)  =  0 


on  6A,        (1.2) 


where  <5A  is  the  boundary  of  the  domain  A;  and  second 


V*U  =  X2U 


subject  to  the  constraint 


d2U 


in  A, 


(1.3) 


U(x,y)  =  Z-Z-   =  0 


on  6A. 


2„2 


(1.4) 


3"n 


In  the  second  category  were  equations  such  as  Poisson's 
equation 


U   +  U   =  f(x,y) 
xx    yy      9a 


subject  to  the  constraint 


in  A, 


(1.5) 


U(x,y)  =  g(x,y) 


and  the  biharmonic  equation 


V*U  =  f(x,y) 


subject  to  the  constraints 


on  6A, 


(1.6) 


in  A, 


(1.7) 


U(x,y)  =  h(x,y) 


on  6A 


(1.8) 


and 


32U 

8n2 


=  q(x,y) 


on  6A. 


(1.9) 


Problems  in  the  first  category  were  numerically  solved 
by  a  relaxation  method,  the  Rayleigh-Ritz  method,  and 
dynamic  programming  combined  with  Stodola's  method. 


Problems  in  the  second  category  were  solved  by  Galerkin's 
method,  the  Rayleigh-Ritz  method  and  dynamic  programming. 
The  domain  used  throughout  this  paper  for  comparisons 
was  the  unit  square,  with  the  constraint 

U(x,y)  =0  on  6A.       (1.10) 

Computations  were  also  carried  out  for  L-shaped  and  trian- 
gular regions  and  for  other  boundary  conditions. 

In  the  relaxation  method  and  dynamic  programming  method 
in  solving  eigenvalue  problems  the  initial  estimates  to 
the  eigenfunctions  had  a  special  property.   The  first  ap- 
proximation to  Ui  was  picked  so  that  it  had  no  negative 
value  in  the  domain.   The  approximation  to  U2  was  picked 
so  that  it  had  negative  and  positive  values  in  the  domain 
and  in  addition  it  was  orthogonal  to  Ui.   The  initial  ap- 
proximations to  U3  was  picked  in  the  same  way  as  the  one 
for  U2  except  that  it  had  to  be  orthogonal  to  both  Ui  and 
U2.  It  may  be  difficult  to  make  a  suitable  choice  for  some 
of  the  higher  modes. 


II.   RELAXATION  METHOD 

For  many  years  relaxation  techniques  have  been  used 
to  solve  differential  equations  with  and  without  the 
aid  of  computers.   They  are  basically  iterative  proce- 
dures in  which  a  new  approximation  is  obtained  from  a 
previous  approximation  and  its  residuals . 

A.   DERIVATION  OF  EQUATIONS 

In  this  section  a  typical  problem  is  posed  and  solved 
by  a  relaxation  method. 

Suppose  the  problem  to  be  solved  has  the  following 
form 

Z,.  =  Z   +  Z  in  A,  (2.1) 

tt    xx    yy  " 

where  the  unknown  function  Z  must  satisfy  the  differential 
equation  in  a  simply  connected  region  A  in  the  xy  plane, 
and  for  t>0.   In  addition  the  function  Z  is  required  to 
vanish  at  points  on  the  boundary  6 A  of  the  region  A, 

Z(x,y,t)  =0  on  6A.  (2.2) 

A  typical  problem  is  that  of  a  vibrating  membrane.   The 
function  Z(x,y,t)  denotes  the  vertical  displacement  of 
the  membrane . 

Now  assume  that  the  displacement  has  a  representation 
of  the  form 

Z(x,y,t)  =  T(t)  U(x,y).  (2.3) 


When  Eq.  (2.3)  is  combined  with  Eq .  (2.1)  the  variables 
may  be  separated  and  a  new  equation  is  obtained  of  the 
form 

mil     P    +  U   ) 

T^=   xx  u  yy  .  _,2 >  (2>4) 

This  is  equivalent  to  two  equations 

T"  +  A2T  =  0  (2.5) 

and 

U   +  u   =  -A2U,  (2.6) 

xx    yy       » 

each  with  a  parameter  A.   Further  U  must  satisfy  the 
boundary  condition 

U(x,y)  =0  on  6A.  (2.7) 

This  is  a  typical  eigenvalue  and  eigenfunction  problem: 

to  find  the  values  X ,  or  A  ,  and  the  associated  functions 

n 

U,  or  U  ,  satisfying  Eq .  (2.6)  and  the  constraint  (2.7). 

The  values  An,  Ai  £  A2  i  A3  <  •••,  are  called  the 

eigenvalues.   With  each  eigenvalue  is  an  eigenfunction 

U  .   In  this  problem  the  eigenfunctions  associated  with 
n  r  ° 

different  eigenvalues  are  orthogonal  on  the  region  A. 

Each  eigenfunction  U  is  also  called  a  mode. 
&  n 

Consider  a  thin  elastic  membrane  of  a  particular  form, 
such  as  a  very  thin  uniform  sheet  of  rubber.   Assume  that 
the  membrane  is  made  fast  at  the  boundary,  while  it  is 
tightly  stretched  over  the  region  with  uniform  tension. 
Also  assume  that  damping  is  negligible.   Then  if  an  interior 
region  of  the  membrane  is  pushed  in  a  direction  perpendicular 


10 


to  the  plane  of  equilibrium,  it  becomes  distorted  into  a 
curved  surface.   The  resulting  area  can  be  computed  as 

S  =  ///l  +  U  2  +  U  2  dxdy,  (2.8) 

A 

Assume  that  U  and  U  are  very  small;   Eq.  (2.8)  becomes 
x      y 

approximately 

S  =  //(I  +  %U  2  +  3§U  2)  dxdy.  (2.9) 

A 

The  increase  in  the  area  of  the  membrane  due  to  the  distor- 
tion is  therefore  approximated  by 

6S  =  JJ(1  +  hV    2  +  ^U  2)  dxdy  -  //  dxdy    (2.10) 
x      y 

A  A 

=  h   //(Ux2  +  Uy2)  dxdy. 
A 

Hence  the  potential  energy  of  the  membrane  in  the  deflec- 
ted position  is 

PE  =  v/2  //(U  2  +  U  2)  dxdy,  (2.11) 

a     y 

A 

where  v  is  the  tension,  assumed  to  be  constant  over  the 
region. 

Now  consider  any  particular  eigenfunction  or  mode. 
It  follows  from  the  solution  of  Eq .  (2.4)  that  the  deflec- 
tion is  a  periodic  function  of  time  and  may  be  expressed 
in  the  form 

Z(x,y,t)  =  U(x,y)  sin  At ,  (2.12) 

except  for  a  phase  shift.   Thus  Eq .  (2.11)  can  be  rewritten 
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as 

PE  =  3g//(U  2    +  U  2)    dxdy      sin2    Xt.  (2.13) 

A 

The  maximum  value  of  the  potential  energy  is 

PE    =  v/2  f f (U  2  +  U  2)  dxdy.  (2.14) 

max       J  J       x  y      J 

A 

The  kinetic  energy  of  an  element  dm  =  pdxdy  of  the  membrane 
is 

1-2pdxdy(Ut)2    =   1-2pdxdy(U2x2    cos2xt),  (2.15) 

where  p  denotes  the  mass  per  unit  area  of  the  membrane. 
Therefore,  the  kinetic  energy  of  the  vibrating  system  is 

KE  =  h\2p    //U2dxdy  cos2Xt,  (2.16) 

A 

and  the  maximum  value  of  the  kinetic  energy  is 

KEmax  =  h\*p    //UMxdy.  (2.17) 

A 

If  it  is  assumed  that  the  energy  is  constant  then,  the 
maximum  values  of  the  potential  and  kinetic  energy  are 
equal  for  individual  modes,  and  therefore 

^X2p  JfU2dxdy  =  v/2  JJ(U2  +  U  2)  dxdy     (2.18) 
*  •  x      y 

A  A 

or, 

v  JJ(Ux2  +  Uy2)  dxdy  (2.19) 

X2  =  —A 

p  //U2dxdy 

A 

12 


The  first  eigenfunction  Ui  is  the  function  which  minimizes 
this  quotient  and  the  first  eigenvalue  Xj2  is  the  corres- 
ponding value  for  the  quotient.   The  next  eigenfunction 
is  the  function  U2  which  minimizes  the  quotient   in  the 
space  of  functions  orthogonal  to  Ui ,  and  X22  is  the  corres- 
ponding value  of  the  quotient.   The  third  eigenfunction 
U3  minimizes  the  quotient  in  the  space  of  functions  ortho- 
gonal to  Ui  and  U2,  etc.   The  set  Ui ,  U2  ,  ...  is  unique 
except  that  if  two  or  more  eigenfunctions  have  the  same 
eigenvalue,  they  may  be  replaced  by  linear  combinations 
of  themselves.   It  is  seen  that  Xi  <  X2  <  X3  ... 

If  the  mode  U  is  known,  Eq.  (2.19)  can  be  used  to  obtain 
X2.   An  interesting  fact  is  that  a  rather  poor  approximation 
to  the  first  mode  Ui ,  chosen  to  satisfy  the  appropriate 
boundary  condition,  will  yield  a  surprisingly  good  estimate 
for  X12.  This  result  is  apparently  due  to  Lord  Rayleigh, 
and  the  quotient  is  sometimes  called  Rayleigh' s  quotient. 

B.   COMPUTATIONAL  ROUTINE 

In  this  section  the  results  of  section  A  will  be  used 
to  develope  a  method  to  solve  Eq .  (2.6). 

If  Eq.  (2.6)  is  expressed  as  a  difference  equation 
then  it  may  be  written  as 

(U4J,  ,  +  U.  -  .  -  2U.  ,)    (U.  ....  +  U.  .  .  -  2U.  .) 
1+1, j     i-l, J     i,j      i,J+l    i,J-l     i9J 

h2  k2 

=  ->2u 

A   i,j,  (2.20) 
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where  h  and  k  denote  the  x  and  y  mesh  size  respectively.   If 
the  mesh  sizes  are  equal  then  Eq.  (2.20)  can  be  rewritten  as 

(U..-  ,  +  U.  ...  +  U.  .  .  +  U.  .  .)     (2.21) 

l+l, j    i,j+l    i-l, J    i,j-l 

u.  .  =  

ljJ  4  -  X2h2 

The  procedure  used  to  solve  the  partial  differential 
equation  was  to  pick  a  function  U°  which  satisfied  the  given 
boundary  conditions  as  a  first  approximation  to  U(x,y). 
This  function  need  not  be  a  very  good  approximation  to  U, 
and  in  fact  step  functions  were  sometimes  used. 

From  Eq .  (2.19)  an  approximation  to  X2  was  obtained 
by  assuming  that  the  approximation  U°  was  the  desired  func- 
tion U.   With  this  first  approximation  to  A2  Eq.(2.21) 
was  used  to  obtain  an  improved  estimate  of  U.   The  form 
of  Eq.  (2.21)  used  to  obtain  the  vth  approximation  was 

v-1      v-1      v        v 

(U. .-  .  +  U.  ...  +  U.  .  .  +  U.  .  .  )    (2.22) 
v      i+l, J     i,J+l    i-l, J     i,J-l 

U.    =  • 


i,J  (4  -  X2h2) 


in  which  the  v-lst  estimate  of  X  was  used.   By  alternating 
Eq.  (2.19)  and  Eq .  (2.22)  a  close  approximation  to  X2  and 
U  were  obtained.   The  solution  converged  to  the  smallest 
eigenvalue  \1    and  the  associated  eigenvector  Ux . 

C.   HIGHER  ORDER  EIGENVALUES 

In  this  section  the  method  is  extended  to  find  the 
larger  eigenvalues  and  eigenfunctions . 

To  obtain  the  larger  eigenvalues  and  eigenvectors  the 
same  procedure  as  that  in  section  B  was  used  except  that 
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additional  equations  were  added  to  force  the  eigenfunctions 
to  be  orthogonal  to  those  already  found.   Define  U.  for 
i  =  lj  2,...j  n  to  be  the  ith_  eigenfunction,  associated 
with  X. .   All  higher  eigenfunctions  must  be  orthogonal  to 
the  lower  ones  obtained.   Orthogonalization  was  accomplished 
by  subtracting  out  multiples  of  the  lower  eigenfunctions 
already  obtained.   This  was  effected  by  expressing  U  as 

U       =U       -  I   CU  (2.23) 

1  new    x  "old   j=l  J  J 

where 


J/U1+1    U.  dxdy 

c>  _  A     °ld ,    j=l,2,...,i.     (2.24) 

J    //U.2  dxdy 
A  J 

For  each  eigenfunction  desired  a  different  initial 
approximating  function  was  used.   The  method  gave  the  eigen- 
values and  associated  eigenfunctions  in  numerically  increas- 
ing order.   The  method  was  very  simple  and  effective  for 
lower  eigenvalues . 

The  fault  of  the  method  is  that  convergence  was  poor 
and  computation  times  were  large  if  the  number  of  intervals 
in  both  directions  was  large. 
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III.   DYNAMIC  PROGRAMMING  [8,10] 

The  method  of  dynamic  programming  has  recently  been 
applied  to  the  solution  of  partial  differential  equations, 
[lOJ .   The  difference  equation  may  be  regarded  as  leading 
to  a  system  of  linear  equations  with  a  large  number  of 
unknowns.   The  method  effectively  reduces  the  number  of 
unknowns  involved  so  that  a  number  of  systems  are  to  be 
solved,  each  one  of     much  lower  dimension.   In  one  case, 
for  example,  using  a  grid  with  N+l  intervals  in  the  region 
in  each  direction,  N  systems  each  with  N  unknowns  are  solved, 
rather  than  one  system  with  N  -  squared  unknowns. 

In  section  A,  the  method  is  applied  to  the  solution  of 
Poisson's  equation  over  a  rectangle,  following  the  paper  of 
Angel  Cl0"3 .   In  section  B,  the  method  is  extended  to  a  re- 
lated eigenvalue  problem  by  combining  it  with  Stodola's 
method.   In  section  C,  the  method  is  extended  to  the  solu- 
tion of  the  biharmonic  equation,  by  applying  dynamic  pro- 
gramming twice.   Finally,  in  section  D,  the  method  is 
applied  to  the  eigenvalue  program  involving  V^U  =  X2U,  again 
by  combining  the  method  with  Stodola's  method. 

A.   DYNAMIC  PROGRAMMING  FOR  PARTIAL  DIFFERENTIAL  EQUATIONS 
Consider  the  solution  of  Poisson's  equation 

U   +  U   =  f(x,y)  in  A,         (3-D 

xx    yy      5J  s 

where  U  =  U(x,y)  is  subject  to  given  boundary  conditions 
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such  as 


U(x,y)  =  g(x,y) 


on  6 A. 


(3.2) 


It  may  be  noted  that  Eq.  (3.1)  is  the  Euler  equation 
associated  with  the  variational  problem 

I(U)  =  Min  //(U  2  +  U  2  +  2fU)  dxdy  (3.3) 

U   A   X      y 

where  the  function  U  is  chosen  from  the  class  of  functions 
with  first  partial  derivatives  belonging  to  L  over  A,  and 
satisfying  the  boundary  conditions  of  (3.2)  on  <Sa. 

Let  the  region  A  be  discretized  by  choosing  n+1  and 
m+1  equally  spaced  points  in  the  independent  variables 
x  and  y  respectively.   Then  Eq .  (3.3)  may  be  rewritten  in 
a  discrete  version  with  equal  intervals  as 


n   m 


I(U)  =  Min    I        I 
U.  .  i=l  j=l 


(U.  .  -  U.  .  J2 


(3.4) 


+  (U.  .  -  U.  .  .)2  +  2f.  .  U,  .h2 
i,J    i-l, J       i,J   i,j 


where  (u    ),  (U    },  {U    },  and  {U.  }   are  determined 
w,j      -L,u     fi » j         i,m 

from  the  boundary  conditions  of  Eq .  (3.2).   Now,  if  all  terms 
in  Eq.  (3.4)  involving  only  boundary  values  were  removed, 
while  not  affecting  the  solution,  a  more  convenient  form 
of  Eq.  (3.4)  is  obtained 


n 
I(U)  =  Min    I 

U.  .  i-l 

i,J 

m-1 


m 


I 


&    CUi,J  "  Ui,J-l} 


(3.5) 


+  I      2h2f.  ,U.  .  +  I      (U.  .  -  U.  ,  .) 
j=l      i.J  i,J    jt1        i,J     i-l, J 


m-1 
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In  vector  notation  Eq.  (3.5)  may  be  rewritten  in  the  form 


n 


I(U)  =  Min  I    (<QUR,UR>  +  <rR,UR>  +  S 


UR   i=l 


R>    R 


R*  R 


R 


(3.6) 


+  <Ur  -  Un_n,  UR  -  UL  ,  >  +  <  fR,UR>) 


'R  -  UR-1>  UR  "  UR-1 


R»  R 


In  this,  UD  =  (U„  , . . .  .  U_    )  .   This  relation  now  defines 

1         m-1 

a  symmetric   matrix  Q,  vectors  rRJ  and  fR,  and  a  scalar 
SRby 


Q  ■  {q,  ,•}.  where  q   . 

-1-  J  J  ±  SO 


2       i  =  J 

-1  |i-j|  =  1 
0  otherwise 


(3.7) 


rR  =  {rR,j}>  Where  rR,j 


1   -2UR,o   J  "  1 

-2UR)m  J  "  "-1 

0  otherwise 


fR  "  {2h2fR,j> 

2  2 

sR  =  UR    +  UD 
R    R,o     R,m 

The  notation  as  in  <rR,  UR> ,  denotes  the  scalar  product  of 
the  two  corresponding  vectors.   The  matrix  Q  remains  constant 
while  rR  and  SR  are  functions  of  the  boundary  conditions 
only. 

Now  in  order  to  solve  Eq .  (3.6)  a  sequence  of  dynamic 
programming  problems  was  considered.   Let 


n 

FR  (V)  =  Min_    I    (<QU. ,  U.>  +  <r, ,  U.> 
UR...Un_1i=R 


(3.8) 


+  s 


.  +  <U.  -  U.  ..  U.  -  U.  ->  +  <f . .  u.>) 
1     1    l-l5   1    1-1      l'   1 
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where  UD  ,  =  V   and  U   is  given  by  Eq.  (3.2).   Now  Eq .  (3.8) 
ti— i  n 

can  be  rewritten  as 

FR(V)  =  Min  (<QURJUR>  +  <FR,  UR>  +  SR         (3-9) 
UR 

+  <ur  _  v,UR  -  V>  +  <fR  ,  UR>  +  <QUR+i,  UR+1> 

+  ...  +  sn  +  <7n  -  trn_13  un  -  vnmml>  +  <Fn,  un>). 

This  may  be  done  since  the  minimization  over  UR+ -,,..., 

U  ,  can  be  commuted  with  the  minimization  over  UD .   This 
n-1  R 

is  a  common  technique  of  dynamic  programming  based  upon  the 
principle  of  optimality,  [8].   This  allows  Eq .  (3-9)  to  be 
rewritten  as 

r 


PD(V)  =  Min 


<QUR,  UR>  +  <rR,  UR>  +  SR        (3.10) 


+  <UR  -  V,  UR  -  V>  +  <fR,  uR>  +  pr+1(ur)  . 

The  final  equation  is 

P  (V)  =  <QU  ,U>  +  <r,U>  +  S  (3.H) 

n        n*  n      n*  n     n 

+  <U   -  V,  U   -  V>  +  <f  ,  u  >, 
n    '   n  n}   n  ' 

and  U   is  known  from  the  boundary  conditions  Eq .  (3.2). 
Since  FR(V)  is  quadratic  in  V  Eq .  (3.11)  may  be  rewritten 
in  the  form 

FR(V)  =  <AR"V,  V>  +  <bRV>  +  CR.  (3.12) 

Now  by  substituting  from  Eq .  (3.12)  into  Eq .  (3.10)  and 
then  differentiating  with  respect  to  UR  an  expression  for 
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UR  is  obtained  of  the  form 


UR  =  (I  +  Q  +  AR+1) 


-1 


V  - 


bR+l  +  rR  +  fR 


.  (3.13) 


This  is  obtained  by  substituting  related  Eq.  (3.13)  into 

Eq.  (3.10)  and  then  combining  Eq .  (3-10)  with  Eq .  (3.12), 

the  various  quantities  in  Eq.  (3.13)  are  defined  by  these 
steps  as 


Ar  =  i  .  (i  +  Q  +  Ar+1) 


-1 


(3.1*0 


bR  =  (I  +  Q  +  AR+1)-X  (bR-1  +  rR  +  fR) 


CR  ~  CR+1  +  SR  " 


(I  +  Q  +  AR+1) 


-1 


<  bR+l  +  rR,  Vl  +  rR  > 
2  2 


(3.15) 
(3.16) 


with  initial  values  determined  from  Eq .  (3.11)  as 


A   =  I, 

n    ' 


(3.17) 


bN  =  "2Un> 


(3.18) 


C   =  <(I  +  Q)  U  ,  U  >  +  <r   -U>+S.       (3.19) 
n  nJn      n    n     n 

The  matrix  (I  +  Q  +  AR  +  ,  )  is  nonsingular,  1  .  Thus  it  has 
inverses  which  may  be  computed  beforehand,  since  it  is  de- 
pendent only  upon  the  type  of  operator. 

Due  to  the  fact  that  only  the  values  of  U^  are  desired 

n 

the  quantities  CD  need  not  be  calculated. 

n 

The  procedure  is  to  calculate  the  quantities  in 
Eq.  (3.14)  repeatedly  until  A2  and  b2  are  obtained.   Then 
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U  is  found  from  Eq.  (3.13)  with  V  =  U  ,  which  is  known 
from  the  boundary  conditions.   Next,  Eq.  (3.13)  is  repeated- 
ly solved  using  the  stored  values  of  AR  and  bR,  and  the 
last  value  of  UR  as  V. 

Thus,  the  problem  was  solved  by  n-1  inversions  of  sym- 
metric matrices  of  order  m-2.   While  these  matrices  may  be 
large  there  are  efficient  computer  routines  available  to 
determine  the  inverses.   Once  the  inverses  are  found,  they 
may  be  stored  for  future  use,  since  they  are  based  only  on 
the  geometry  of  the  region  A.   Thus  for  several  problems 
over  the  same  region  the  inverses  may  be  entered  into  the 
program  as  data.   Also  they  have  the  property  that  if  less 
than  n+1  grid  points  are  required  a  reduced  number  of  the 
matrices  may  be  used. 

B.   DYNAMIC  PROGRAMMING  FOR  EIGENVALUE  PROBLEMS 

In  the  first  section  of  this  chapter  it  was  seen  that 
dynamic  programming  could  be  used  to  solve  partial  differ- 
ential equations.   In  this  section,  the  program  is  modified 
to  solve  a  related  eigenvalue  problem  by  Stodola's  method. 
Consider  the  problem  of  the  vibrating  membrane  consi- 
dered in  Chapter  II.   The  differential  equation  for  the 
eigenvalues  and  eigenfunctions  is 

U    +  U    =  -A2U  in  A,  (3-20) 

xx    yy 

where  U  =  U(x,y)  is  subject  to  the  boundary  condition 

U(x,y)  =0  on  6A.         (3.21) 
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This  problem  was  solved  by  two  related  methods.   In 
both  Eq.  (3.20) is  regarded  as  a  special  case  of  Eq.  (3.1) 
in  which 


and 


f(x,y)  =  -A2U(x,y)  in  A  (3.22) 

f(x,y)  =0  on  6A         (3.23) 

1.   First  Routine 

The  difficulty  here  was  that  the  function  f  was 
known  only  on  the  boundary  and  \   was  unknown.   The  first 
step,  to  overcome  this  obstacle,  was  to  choose  a  function 
u°,(x,y) which  satisfied  the  given  boundary  conditions.   Then 
an  estimate  of  A2  was  obtained,  say  (a  )2  ,  by  using  Ray- 
leigh's  formula.   Next  a  new  estimate  of  U,  say  \J1(x)y),   was 
obtained  using  dynamic  programming  to  solve  the  equation 

Uxx  +  Uyy  =  -(*°)2  U°  =  f°(x,y).  (3-2*) 

By  repeating  this  sequence  of  Rayleigh's  formula  and  dyna- 
mic programming,  a  good  approximation  to  the  minimum  eigen- 
value and  the  associated  eigenfunction  was  obtained.   The 
next  two  eigenvalues  and  vectors  were  also  obtained  by  the 
process  of  forcing  the  higher  eigenfunctions  to  be  orthogo- 
nal to  ones  already  obtained  as  was  done  in  relaxation. 
The  inverse  matrices  used  in  the  routine  were  calculated 
in  determining  the  first  eigenvalue  and  eigenfunction; 
they  were  then  stored  so  that  it  was  not  necessary  to  recal- 
culate them  for  subsequent  eigenvalues. 
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2.   Stodola's  Method 


The  second  form  is  more  like  the  usual  form  of  Sto- 
dola's method  and  hence  a  brief  of  Stodola's  technique  is 
given  first . 

An  initial  function  V0  satisfying  the  boundary  con- 
ditions is  selected.   It  may  be  considered  to  be  of  the 
following  form 


Vo  =  aiUi  +  a2U2  +. .  . 


(3.25) 


where  ai  is  not  equal  to  zero.   For  convenience  V  was 


normalized;   the  L^  norm  of  V0, 

MAX  |V0(x,y)|  =  | |V0| |; 

A 


V, 


,  is  defined  as 


(3.26) 


and   |V0|   is  set  equal  to  one.   Now  consider 


L  V0  -  ^4  Ui  -  ^-r,  U2 


(3.27) 


and 


L-mV0  = 


-1 


l*~^J 


m 


aiUi  +  a: 


f,    i  2m 

A  i 
A2 


U2  +. . . 


(3.28) 


In  Eq .  (3.28)  it  can  be  seen  that  the  relative  size  of  the 
components  U2  ,  U3,...  are  decreasing  by  a  factor  (Ai/A2)2, 
(\i/A3)2,...  respectively.   After  a  few  applications  of  the 
operator,  L   ,  to  VQ,  the  leading  term  will  dominate. 
Thus  the  functions  obtained  will  approximate  V x  ,  except 
for  a  constant  factor,  and  the  ratio  of  successive  iterates 
will  approach  a  constant,  -1/x2. 
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3 .   Second  Routine 

This  was  applied  as  follows .   Let  Z   be  defined  by 


the  equation  LZ   =  V   . .  so  t 

m    m-1 ' 


hat 


Zm  "  L"Vr  (3-29) 

This  equation  was  solved  by  dynamic  programming.   It  was 

found  convenient  to  normalize  each  iteration.  Let 

-Z 

V  =  ^—   .  (3.30) 

||Z  || 
i  i  mi  i 

Then  the  approximation  can  be  made 

Ai2  =  —i .  (3.3D 

|  |  Z  || 

The  sequence  of  functions  V0 ,  Vi,...  converges  to  Ui . 

The  process  was  terminated  when  successive  approx- 
imations were  sufficiently  close  together,  say, 

I |Vm(x,y)  -  Vm-1(x-y)| |<e,  (3-32) 

where  e  is  a  preassigned  small  positive  number.   The  result- 
ing function  V  is  an  approximation  to  Ui ,   and  Xi  is  ob- 
tained from  Eq.  (3. 3D* 

C.   DYNAMIC  PROGRAMMING  FOR  A  HIGHER-ORDER  OPERATOR 

In  section  A  and  B  of  this  chapter  it  was  shown  how 
dynamic  programming  could  be  used  to  solve  Poisson's  equa- 
tion and  the  vibrating  membrane  program,  respectively. 
In  this  section,  by  another  modification  to  the  routine, 
the  biharmonic  equation  can  be  solved. 


2k 


Assume  the  problem  to  be  solved  is  of  the  form 


V(U)  =  f(x,y) 


in  A, 


(3.33) 


where  U  =  U(x,y)  is  subject  to  the  constraints 

'  U  =  g(x,y)  on  6A        (3.34) 

■ 

32U/3n2  =  P(x,y)  on  6A. 

Clearly  Eq.  ( 3 • 33 )  may  be  rewritten  in  the  form 

V2  *(x,y)  =  f(x,y),  (3.35) 

where 

V2  U(x,y)  =  *(x,y).  (3-36) 

Now,  since  U  and  92U/3n2  are  both  known  as  the  boundary, 
*(x,y)  may  be  approximated  on  the  boundary.   On  a  rectangle, 
for  example,  the  following  relations  determine  $  on  the 
boundary 

*(0,y)  =  -P(0,y)  +  Uyy  (0,y)  (3-37) 

$(n,y)  =  P(n,y)   +  n^    (n,y) 

4(x,0)  =  -P(x,0)  +  U    (x,0) 

A  A 


k  $(x,m)  =  P(x,m)   +  U    (x,m) 

A.  .A. 


The  usual  finite  difference  scheme  may  be  used  to  approxi- 
mate U   and  U   and  thus  the  relations  of  Eq.  (3.37)  may 
xx      yy  ^  * 

be  approximated  by 


$  =    -P  + 

0,j  F0,j 


(Un    , ,-    -    2Un     .    +   Un     .    -  ) 

0,j+l  0,j  9juL±1      (3   38) 

h2 
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=  p  .  + 

n,j    n,j 


(U   . ,.  -  2U   ,  +  U   ,  .  ) 

...1+1     n,J     nj-r    (3.38) 


•-  n  -  -P.  n  + 

1,0      1,0 


(Ui+1.0 


h1 


-  2U.  n  +  U.  .  n) 
1 ,0 l-l  ,0 


h' 


=  P 


(U.j..    -  2U,    +  U.  .   ) 

l+l,  m 1  ,m l-l  ,m 


for  i  =  1,...,  n-1,  j  =  1,...,  m-1, 

and  <±>n  «.  $   n ,  $n   ,  $    are  known  from  the  boundary  con- 
0,0'   n,0 J   0 ,m'   n,m  J 

ditions.   Thus  $(x,y)  is  now  approximated  on  the  boundary. 
First  Eq.  (3. 35)  is  solved  by  dynamic  programming  to  obtain 
the  function  $  in  the  region.   Then  Eq.  (3.36)  is  solved 
by  dynamic  programming  to  obtain  the  desired  solution. 


D.   DYNAMIC  PROGRAMMING  FOR  HIGHER  ORDER  EIGENVALUE  PROBLEMS 

This  method  may  be  extended  to  the  corresponding  eigen- 
value and  eigenfunction  problem  much  as  before.   One  such 
physical  problem  is  that  of  a  vibrating  uniform  plate  with 
hinged  edges. 

Assume  that  the  differential  equation  has  the  form 


V^U  =  X2U, 


and  the  boundary  conditions  are 


u  -  ilE.  0 

8n2 


in  A, 


(3.39) 


on  6A 


(3.40) 


The  technique  of  sections  B  and  C  may  be  applied  in  two 
steps  to  the  solution  of  the  problem.   Let  $  be  defined  as 


$  =  V2U. 


(3.^1) 
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Then  the  two  equations  to  be  solved  are 

V2$n  =  x2Un  (3.42) 

and 

V2un+1  =  $n,  (3.43) 

subject  to  the  boundary  conditions,  Eq.  (3.40). 

In  order  to  start  the  routine  an  initial  estimate  for 
the  function  U(x,y),  say  {U..},  is  made.   The  value  of  X 
is  estimated  as  was  done  in  section  B.   Then  Eq.  (3.^2) 
and  Eq.  (3.^3)  are  solved  by  dynamic  programming  to  get  the 
next  approximation  {U.  .}.   The  same  criterion  for  stopping 
is  used  as  that  in  section  B. 
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IV.   RAYLEIGH-RITZ  METHOD  IX  5, 6, 9] 

The  Rayleigh-Ritz  method  has  been  used  for  many  years 
to  obtain  approximations  to  the  solution  of  partial  differ- 
ential equations  and  eigenfunction  problems.   In  this  method 
the  problem  is  posed  as  a  minimization  problem,  say  invol- 
ving an  integral.   Then  some  linearly  independent  functions 
which  satisfy  the  boundary  conditions  are  chosen.   The  solu- 
tion is  approximated  by  a  linear  combination  of  these. 
Finally  the  coefficients  in  the  approximation  are  chosen 
so  as  to  effect  minimization.   This  leads  to  an  eigenvalue 
problem  involving  symmetric  matrices. 

The  functions  chosen  may  be,  for  example,  polynomials 
of  low  degree,  or  trigonometric  functions.   Assume  that 
homogeneous  boundary  conditions  are  given.   Let  $,  =  $,  (x,y) 
be  n  functions  which  satisfy  these  and  approximate  the  solu- 
tion U  by 

U(x,y)  =  I      C  $  (4.1) 

k=l   K  K 

A.   PARTIAL  DIFFERENTIAL  EQUATIONS 
In  this  section  the  solution  of 

U   +  U   =  f(x,y)  in  A,         (1.3) 

xx     yy       ,l/  J 

is  again  considered.   It  is  the  Euler  equation  associated 
with  minimizing  the  integral 
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//(U  2  +  U  2  +  2fU)  dxdy.  (4.2) 

A   X      y 

Substituting  from  Eq.  (4.1)  into  Eq .  (4.2)  and  integration 
results  in  a  function  which  may  be  written 

I  =  I(Cl,...,  Cn),  (4.3) 

for  functions  of  the  form  (4.1).   The  minimization  of  this 
function  and  an  approximate  solution  to  Eq.  (1.3)  is  thus 
obtained  by  solving 

|I-  =  0,     k  =  1,  2,...,  n.  (4.4) 

The  effectiveness  of  the  procedure  of  course  depends  upon 
the  choice  of  the  approximating  functions  $,  (x,y). 

B.   EIGENVALUE  PROBLEMS 

The  method  is  also  applicable  to  eigenvalue  problems. 
Consider  again  the  problem  in  Chapter  II  of  finding  eigen- 
values and  eigenfunctions  for  the  equation 

U   +  U   =  -A2U  in  A         (2.6) 

xx    yy 

subject  to  the  boundary  condition 

U(x,y)  =0  on  «A.       (2.7) 

2 

In  many  problems  Xi,  the  lowest  eigenvalue,  is  the  minimum 
of  the  ratio  of  two  integrals.   This  fact  was  shown  in 
Chapter  II  and  Eq.  (2.19). 

If  an  approximation  to  U  is  chosen  as  in  section  A,  since 
each  function  satisfies  the  linear  homogeneous  boundary 
condition,  the  sum  shown  in  Eq.  (4.1)  satisfies  it.   If 
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this  sum  is  substituted  into  Eq.  (2.19)  the  resulting  values 

define  an  upper  bound  for  A,  for  all  choices  of  the  C's. 

Further  also,  the  C's  may  be  chosen  to  give  a  least  upper 

bound  over  the  subspace  spanned  by  Si,  $2,...,  $  . 

If  that  approximation  for  U  is  substituted  into  Eq.(2.19) 

the  numerator  and  denominator  become  quadratic  forms  in 

(C  .....  C  )  =  C.   The  numerator  has  the  form 
1  '      n 

I        I      a.        C  C   =  C^A  C  (4.5) 

1=1  j=l   1J  ±   J 


where 


=  //(*   <D   +  •   •   )  dxdy,  (4.6) 

J    A    x  Jx     y  Jy 


and  the  denominator  has  the  form 


n   n 


ttW 


I         I      b.  .C.C.  =  CTB  C  (4.7) 

1=1  j=l   1J  X  J 


where 


b   =  //•   •  dxdy.  (4.8) 

J- J    A  x   J 

These  define  the  matrices  A  and  B.   The  first  eigenvalue 
Xi  minimizes  the  quotient  in  Eq.  (2.19)  and  hence  the  mini- 
mum value  of 

C^A  C/C^B  C 

defines  the  minimum  value  of  the  quotient  in  the  subspace 
spanned  by  $!,,..,  $  . 
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To  find  this  is  equivalent  to  minimizing  the  quadratic 
form 

Xi2  =  Min  C^A  C  (4.10) 

C 

subject  to  the  constraint  that  the  second  quadratic  form 
assumes  the  value  one 

0%  C  =  1.  (4.11) 

The  problem  may  be  solved  in  two  steps.   First  find  the  eigen- 
values and  eigenvectors  of  B.   Let  the  eigenvalues  of  B 
be  a.2,  i=l,  2,...,  n,  and  the  associated  normalized  eigen- 
vectors V.  .   Let  T  be  the  matrix 

i 

T  =  (Vi,  V2,...,  Vn)  diag(  l/a)i  ,.  ..,  1/%).   (4.12) 

Then  the  transformation 

C  =  T  D  (4.13) 

reduces  the  constraint  (4.11)  to  the  form 

C^B  C  =  ^D  =  1.  (4.14) 

condition  (4.10)  becomes 

Xi2  =  Min  5*1  D,  (4.15) 

D 


where 


E  =  T*A  T,  (4.16) 


subject  to  the  constraint  (4.14). 

Hence  the  value  of  Xi2  is  the  smallest  eigenvalue  of 
the  matrix  E. 
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Let  Di  be  the  associated  normalized  eigenvector  of  E. 
Then  the  corresponding  minimizing  function  has  coefficients 
determined  from  Eq.  (4.13) 

Ci  =  T  Di.  (4.17) 

This  value  of  Ai2is  an  upper  bound  for  the  first  eigen- 
value, the  least  upper  bound  in  the  space  of  functions 
spanned  by  $1  ,...,$  .   In  a  similar  way  the  second  eigen- 
value of  E  furnishes  an  upper  bound  to  the  second  eigenfunc- 
tion,  etc. 

For  simple  problems,  at  least,  it  seems  easy  to  choose 
the  $'s  so  that  a  good  bound  for  the  lowest  eigenvalue  re- 
sults.  It  is  not  clear  how  to  make  good  choices  to  get 
good  estimates  for  the  higher  eigenvalues. 
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V.   GALERKIN'S  METHOD  [4,5, 6] 


Sometimes  the  solution  to  boundary-value  problems  in- 
volving partial  differential  equations  can  be  obtained  by 
forming  an  associated  functional  in  the  form  of  a  definite 
integral  which  is  to  be  made  stationary.   For  a  problem  of 
this  type,  the  Rayleigh-Ritz  method  is  often  an  effective 
procedure  for  the  determination  of  an  approximate  solution, 
However,  in  many  instances  it  is  difficult  to  find  this 
functional.   In  such  cases,  the  Galerkin  method  is  often 
effective . 

A.   PARTIAL  DIFFERENTIAL  EQUATIONS 

Consider  the  linear  homogeneous  boundary  value  problem 
of  the  form 

LU(x,y)  =  f(x,y),  (5.1) 

subject  to  linear  homogeneous  boundary  conditions.   The 
symbol  L  stands  for  a  linear  differential  operator,  such 
as  V2  • 

Suppose  that  an  approximate  solution  is  taken  in  the 
form 

Un(x,y)  =  I      Ck$k(x,y),  (5.2) 

k=l 

where  the  coefficients  C, , . . . ,  C   are  constants,  to  be  de- 
termined, and,  as  in  the  previous  chapter,  the  functions 
$,,...,  $   are  picked  to  satisfy  the  homogeneous  boundary 
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conditions.   The  coefficients  are  dependent  upon  the  number 
of  functions  picked  and  therefore,  must  be  recomputed  if 
a  larger  number  of  functions  are  later  chosen.   While  the 
function  U   satisfies  the  boundary  conditions,  it  will  not, 
in  general,  satisfy  Eq.  (5.1).   Thus  there  is  a  residual, 

LUn(x,y)  -  f(x,y)  =  Rn(x,y),  (5.3) 

which  can  be  viewed  as  an  error  or  penalty  function.   Assume 
the  solution  for  U(x,y)  can  be  expressed  by  an  infinite 
complete  series  of  these  linearly  independent  functions 
in  the  form 


U(x,y)  =  I      o   <b.   (x,y).  (5.4) 

k=l  K  K 

Now  U  (x,y)  in  Eq.  (5.2)  represents  a  sequence  of  partial 
sums  which  approximate  U(x,y).   Now  if  the  condition  is 
imposed  that  L(U  )  -  f  be  orthogonal  to  each  function 
$.(x,y)  on  the  demain  A,  the  following  set  of  equations 
are  obtained 


(5.5) 


/J(L(U  )  -  f)    $,  dxdy  =  0,  for  k  =  1,  2,...,  n 
A     n        k 

or  by  use  of  Eq.  (5-3)  (5.6) 

//R  (x,y)  $  (x,y)  dxdy  =  0,  for  k  =  1,  2,...,  n. 
A  n       k 

If  Eq.  (5.6)  is  to  hold  as  n  ■*■  °°  it  follows  that 

lim  Rn  =  0  (5.7) 
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Let  R(x,y)  be  an  arbitrary  function  satisfying  the 
homogeneous  boundary  conditions.   Since  {$,  )  for  k  =  1,  2,... 
forms  a  complete  set  of  functions,  constants  a.,,...,  a  , ... 

can  be  found  such  that 


R(x,y)  =   I   a.$v(x,y). 


k=l 


k  k 


(5.8) 


Thus, 


//   lim  R    (x,y)    n(x,y)    dxdy 
A    n+°° 


=  0, 


(5.9) 


for  any  arbitrary  n(x,y),  so  that  by  the  fundamental  lemma 
of  variational  calculus  Eq.  (5.7)  is  true  "almost  every- 
where".  Suppose  further  that  L(U  )  ■*  L(U);  then  Eq.  (5.1) 
holds,  by  the  use  of  Eq.  (5.7). 

Galerkin's  method  requires  that  the  error  function 
R  (x,y)  be  orthogonal  to  each  of  the  functions  $,  ,  or  that 
Eq.  (5.5)  and  Eq.  (5.6)  hold.   Now  by  substituting  Eq .  (5.2) 
into  Eq.  (5.5)  an  integral  is  obtained  of  the  form 


// 

A 


n 


L(  I  CVS>)  -  f 


k=l 


k  k' 


»1  dxdy  =  0, 


(5.10) 


i  =  1,  2,...,  n. 


This  results  in  a  system  of  n  linear  algebraic  equations 

in  n  unknowns  Cn  ,  .  . .  ,  C  .   Furthermore,  the  system  is  in- 

1 '    '   n  ' 

homogeneous  unless  the  function  f(x,y)  is  orthogonal  to 
each  $.  (x,y).   Eq.  (5.10)  may  be  rewritten 
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I   C.   //L(*  )  *   dxdy  =   //f«  dxdy,      (5.11) 
j=l   J   A    3        K         A   K 

k  "■  _i_  .  (—»...»  r  i  • 

In  matrix  notation  this  becomes 

B  C  =  F,  (5.12) 


where 


C  =  (C1   C2  ...Cn)T  (5.13) 

B  =  {b±J}   i  =  1,  2,...,  n,                (5.14) 

j  =  1,  2, . . . ,  n 

F  =  (f1,...,  fn)T  (5.15) 


and 


b. .  =  J/L($. )  ♦.  dxdy  (5.16) 

!<]      A 

f,  =  //f*  dxdy  (5.17) 

x  A   x 

The  constants   c.  are  obtained  by  solving  the  system  Eq. 

Eq.  (5.12). 

Collocation,   a  convenient  variation.   One  way  to  get 

an  approximate  solution  for  the  C's  is  the  following. 

Choose  n  points  of  the  region  A.   Evaluate  the  terms  in 

the  integral  of  Eq.  (5.10)  at  these  points.   This  yields  n 

equations  for  the  unknowns  C  .   This  method,  called  colloca- 

k 

tion,  is  not  generally  so  accurate  but  it  is  quicker  than 
carrying  out  the  integrations  to  define  the  coefficients 
Eq.  (5.16)  and  Eq.  (5.17) . 


36 


Galerkins  method  is  also  useful  to  obtain  a  direct 
solution  of  variational  problems.   Assume  it  is  desired  to 
minimize  a  functional  of  the  form 

KV)  =  //P(x,y  V,V.V  )  dxdy.  (5.18) 

A         x   y 

It  is  desired  to  find  an  extreme  value  for  the  functional 
subject  to  the  condition  that  V(x,y)  is  prescribed  on  the 
boundary  <5A  of  the  domain  A.   It  is  known  that  if  the  exist- 
ence of  an  extremizing  function  U(x,y)  is  assumed  and  that 
the  function  F  possesses  continuous  derivatives  of  the 
second  order  with  respect  to  its  arguements,  there  results 
the  condition 


JJn(x,y) 

A 


Fu  "  T~  Fu  "  ~  Fu 


dxdy  =  0    (5.19) 


x    x     y 

for  an  arbitrary  function  n(x,y)  which  has  piecewise  con- 
tinuous derivatives  in  A  and  which  vanishes  on  6A.   This 
is  derived  by  considering  a  function  of  the  form 

V(x,y)  =  U(x,y)  +  en(x,y),  (5.20) 

and  differentiating  with  respect  to  e.   If  n(x,y)  is  other- 
wise arbitrary,  then  one  form  of  the  fundamental  lemma  of 
variational  calculus  requires  that 

Fu"FFu  -  f-  Fu  =  °-  C5-21) 

x    x   °y    y 

Suppose  now  that  n(x,y)  is  the  kth  function  <f>,(x,y)  and 
that  an  orthogonality  requirement  is  imposed  as 


37 


'/VFu  -  h  Fux  -  S7  FV  dxdy  =  "•  (5-22) 

for  k  =  1,  2,...,  n. 

Finally,  assume  that  the  function  U(x,y)  which  effects  the 
minimization  can  be  represented  satisfactorily  by  a  finite 
series 

U(x,y)  =  I      c  * ,(x,y).  (5.23) 

n        k=l  K  K 

Eq.  (5.22)  then  defines  a  system  of  n  algebraic  equation 

to  be  solved  for  the  n  unknowns  c,  . . . .  c  .   Thus  Galerkin's 

1      n 

method  is  applicable  in  the  solution  of  variational  problems. 
However,  much  of  its  value  lies  in  the  fact  that  it  is  not 
necessarily  connected  with  a  variational  procedure. 

B.   EIGENVALUE  PROBLEMS 

Suppose  that  it  is  desired  to  solve  an  eigenvalue  prob- 
lem of  the  form 

LU  =  AU(x,y)  in  A,        (5.24) 

and 

U  =  0  on  6A.        (5.25) 

Assume  as  in  the  first  section  that  the  function  U(x,y) 
can  be  approximated  in  the  form  of  Eq .  (5.23).   The  problem 
is  solved  by  considering  equations  of  the  form 


'/ 


L(U  )  -  AU 
ny     n 


$.  dxdy  =  0  (5.26) 

j 


for  j  =  1 ,  2 , . . . ,  n . 
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Now  by  substituting  for  U  the  functions  $.  a  system  is  ob- 

n  AC 

tained  of  the  form 


I   (a   -  XY.   )C   =  0  (5.27) 

j=l   1,J     X,J   J 

for  i  =  1,  2,...,  n, 


where 


a.    =  //L($.)  •  dxdy  (5.28) 

Y.  ,  =  //«-*,  dxdy  (5.29) 

A 

This  has  nontrivial  solutions  for  the  C's  if,  and  only  if, 

A  =  | a,   -  Ay,   |  =  0,  (5.30) 

where  A  is  the  determinant  of  the  matrix.   Thus  the  values 
of  X.  may  be  found  by  solving  the  characteristic  equation 
(5.30).   For  each  eigenvalue  X.  there  corresponds  a  system 
of  equations  (5.28)  to  be  solved  for  the  eigenvector 

X  K,j  "  *iYk>j>  Ck  -  0  (5.31) 

for  j  =  1 ,  2  , .  .  .  ,  n . 

t 

Now  for  each  eigenvalue  \.  this  system  of  n  homogeneous 
equations  may  be  solved  to  give  the  values  C,  ,  where 
this  coefficient  corresponding  to  the  ith  eigenvalue. 

Thus  the  eigenvalues  are  obtained,  with  their  correspond- 
ing eigenvectors.   The  functions  $.  should  be  chosen  to 
reflect  whatever  characteristics  the  solution  is  felt  to 

have. 
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VI.   DISCUSSION  OF  VARIOUS  METHODS 


FOR  SOLVING  EIGENVALUE  PROBLEMS 


In  Chapters  II,  III  and  IV,  three  methods  were  developed 
for  finding  the  first  three  eigenvalues  and  eigenfunctions . 
All  three  gave  satisfactory  values  for  the  eigenvalues  and 
eigenfunctions,  and  satisfactory  times  for  the  rather  simple 
problems  considered  here.   These  methods  were  used  to  obtain 
solutions  of  the  following  two  problems: 

in  A,        (6.1) 


U   +  U   =  -  X2U 
xx    yy 


subject  to  the  constraint 


and 


U(x,y)  =  0 


V*U  =  X2U 


subject  to  the  constraints 


U(x,y)  = 


32U 

8n2 


=  0 


on  6A; 


in  A, 


on  6A. 


(6.2) 


(6.3) 


(6.4) 


In  this  chapter  the  computation  and  numerical  results  are 
discussed  and  compared. 

In  all  of  the  methods  a  set  of  functions  was  needed. 
In  the  Rayleigh-Ritz  method  these  were  the  basis  for  the 
approximating  functions;   in  the  other  two  methods  they 
were  the  initial  estimates  of  the  functions.   The  ones  usu- 
ally chosen  were 
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Ui°  =  (x  -  x2)(y  -  y2),  (6.5) 

U2°  =  (x  -  x2)(y  -  y2)0-2  -  x),  (6.6) 

U3°  =  (x  -  x2)(y  -  y2)(^  -  y).  '  (6.7) 

Step  functions  were  also  used  as  first  estimates  in  the 
iterative  methods;   these  increased  the  number  of  iterations 
required  some,  but  not  much,  particularly  for  the  higher 
modes . 

Usually  the  range  of  the  independent  variable  was  di- 
vided into  ten  equal  sub-intervals,  in  going  to  a  differ- 
ence equation.  In  the  Rayleigh-Ritz  method  this  did  not 
yield  sufficient  accuracy  and  it  was  found  necessary  to  go 
to  forty  intervals.  Twenty-five  intervals  were  also  used 
in  some  computations;  the  intermediate  number  was  chosen 
because  of  storage  requirements  in  dynamic  programming. 

In  the  iterative  procedures  some  convergence  or  stop 
criterion  was  needed.  When  a  relaxation  procedure  was  used 
together  with  Rayleigh's  formula  for  estimating  the  eigen- 
value, computation  was  terminated  whenever  the  eigenvalue 
did  not  decrease  by  at  least  0.002.   In  Stodola's  method, 
the  routine  was  terminated  whenever  the  norm  of  the  change 
in  the  function  U(x,y)  was  less  than  0.002. 

The  relaxation  method  required  the  least  time  for  this 
simple  problem.   Most  of  the  time  required  for  dynamic  pro- 
gramming was  spent  in  inverting  the  matrices.   Dynamic  pro- 
gramming yielded  a  very  accurate  approximation  to  the 
eigenfunction. 
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For  some  reason  it  was  necessary  to  use  forty  intervals 
to  get  the  desired  accuracy  with  the  Rayleigh-Ritz  method. 
When  ten  intervals  were  used,  the  eigenvalues  of  the  matrix 
E  were  significantly  too  small,  apparently  due  to  errors 
in  the  integration  and  transformation  routines.   There  were 
at  least  two  other  disadvantages  of  the  Rayleigh-Ritz  method 
First  it  was  tedious  to  program.   Second,  there  may  be  some 
difficulties  in  choosing  the  functions  $.,  particularly 
if  some  of  the  higher  modes  are  desired.   The  results  of  a 
set  of  computations  are  given  in  Computer  Output  5  ,  in- 
cluding the  three  eigenvalues,  the  associated  coefficients 
and  the  values  of  the  corresponding  functions  at  various 
points.   The  routine  is  shown  in  Computer  Program  5«   The 
necessary  matrix  transformations  and  solutions  for  the 
eigenvalues  were  carried  out  using  programs  TRED2  and  TGL2 
respectively,  [73. 

The  simple  relaxation  method  of  Chapter  II  had  the  ad- 
vantage of  being  the  simplest  to  program  and  to  run.   It 
had  the  disadvantages  that  terminal  convergence  was  slower 
than  in  the  method  of  dynamic  programming  and  if  a  large 
number  of  points  were  involved  the  computing  time  increased 
greatly.   The  results  of  a  set  of  computations  are  given 
in  Computer  Output  1.   The  routine  is  shown  in  Computer 
Program  1. 

Dynamic  programming  had  the  advantage  of  yielding  very 
accurate  values  in  a  small  number  of  iterations.  The  dis- 
advantages were  that  it  was  relatively  difficult  to  program 
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and  that  it  took  quite  a  bit  of  time  to  generate  the  inver- 
ses of  the  matrices  required.   Computer  Output  7   shows  the 
results  by  this  method  and  the  program  used  is  in  Computer 
Program  6 . 

The  values  of  the  eigenvalues  obtained  by  the  three 
methods  are  compared  in  Table  I.   In  Table  I  are  the  eigen- 
values for  the  first  differential  equation,   the  difference 
equation  with  ten  intervals,  together  with  the  results  of 
the  computations,  the  number  in  parenthesis  by  an  entry 
indicates  the  number  of  intervals  used  in  the  computation. 
The  eigenvalues  for  forty  intervals  are  intermediate  between 
those   for  ten  and  those  for  the  differential  equation. 


Differential 

Equation 

Difference 
Equation (10) 

Rayleigh 

Relaxation 

Dynamic 
Prograrrming 

4.44289 

4.42463 

4.4110    (25) 

4.42486(10) 

4.42442(10) 

4.44751(40) 

4.44611(40) 

4.43753(25) 

7.02482 

6.92714 

7.23029(40) 

6.92736(10) 

6.92717(10) 

7.02482 

6.92714 

7.24707(40) 

6.92736(10) 

6.92717(10) 

Table  I. 
Comparison  of  First  Three  Eigenvalues  for  Eq.  (1.1) 

The  three  methods  gave  close  agreement  for  the  first 
eigenvalue  but  tended  to  differ  on  the  next  two. 

For  the  differential  equation  of  higher  order,  only 
relaxation  and  dynamic  programming  were  compared;   these 
comparisons  are  shown  in  Table  II  using  ten  intervals. 
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Differential 
Equation 

Difference 
Equation(lO) 

Relaxation 

Dynamic 

Programming 

19.739227 

19.577361 

19.60767 

19.57777 

49.348040 

47.985220 

48.00555 

47.98473 

49 .  3^8040 

47.985220 

48.02386 

47.98474 

Table  II. 
Comparison  of  First  Three  Eigenvalues  for  Eq.  (1.3) 

Computing  times  were  similar,  around  twelve  seconds 
for  each.   Computer  Output  2  shows  the  results  for  the  re- 
location method,  and  the  program  is  Computer  Program  2. 
Computer  Output  9  shows  the  results  for  dynamic  programming 
and  the  program  is  Computer  Program  6. 
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VII.   DISCUSSION  OF  THE  SOLUTION  OF  A 
PARTIAL  DIFFERENTIAL  EQUATION  BY  VARIOUS  METHODS 


In  Chapters  III,  IV,  and  V,  three  methods  were  developed 
for  solving  partial  differential  equations.   These  methods 
were  used  to  obtain  solutions  of  the  following  problems: 

Uxx  +  Uyy  =  2(x*  +  y2  -  x  "  y)     in  A>     (7.1) 


subject  to  the  constraint 


U(x,y)  =  0 


and 


V"U  =  8 
subject  to  the  constraints 
U(x,y)  =  0 
82U 


3n: 


=  2(y  -y2) 


on  6A;    (7.2) 


in  A,     (7.3) 


on  6A 


for  x  =  0 


(7.4) 


82U 
8n2 


=  -2(y  -y2) 


for  x  =  1 


82U 


=  2(x  -x2) 


for  y  =  0 


32U 
8n2 


=  -2(x  -x2) 


for  y  =  0 


In  this  chapter  the  computation  and  numerical  results  are 
discussed. 
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In  Galerkin's  method,  an  approximation  for  the  function 
U(x,y)  satisfying  the  constraint  Eq.  (7.2)  was  made  by  choos- 
ing suitable  functions  $i(x,y),  $2(x,y)  and  *3(x,y)  to  satis- 
fy the  constraint  in  Eq.  (7.2).   The  approximation  for  U(x,y) 
was 

U(x,y)  =  Ci*i(x,y)  +  C2$2(x,y)  +  C3$3(x,y)    (7.5) 

=  (x  -x2)(y  -y2)(Ci  +yxC2  +x2y2C3). 

The  values  of  Ci,  C2  and  C3  were  obtained  by  techniques 
described  in  Chapter  V.   Computer  Output  3  shows  the  values 
of  Ci,  C2,  C3  and  the  function  U(x,y)  at  various  points. 
In  this  method  the  region  A  was  sub-divided  into  ten  equal 
sub-intervals  for  each  independent  variable  for  the  inte- 
gration routine.   A  second  approximation  for  U(x,y)  was  made 
for  this  method  as 

U(x,y)  =  Ci*i  +  C2<S>2  +  C3$3  (7.6) 

=  (x  -x2)(y  -y2)(xC!  +  yC2  +  XyC3). 

The  purpose  of  this  was  as  follows.   There  generally  is  some 
skill  and  art  involved  in  choosing  the  functions  $-,,...,  $ 
well.   In  fact  in  the  first  choice  $-  is  actually  the  desired 
function.   The  second  set  of  functions  was  chosen  so  as  to 
get  some  feel  for  the  consequences  of  a  poor  choice  of  the 
$.'s.   The  values  of  Clt    C2,  C3  and  the  function  at  various 
points  is  shown  in  Computer  Output  10. 

The  numerical  solution  of  Eq .  (7.1)  by  dynamic  program- 
ming with  the  constraint  of  (7.2)  yielded  the  values  in 
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Computer  Output  6.   Again  the  region  was  divided  as  was 
done  in  Galerkin's  method. 

In  the  Rayleigh-Ritz  method,  an  approximation  for  the 
solution,  U(x,y),  satisfying  Eq.  (7.1)  was  made  by  choosing 
suitable  functions  $i(x,y),  $2(x,y)  and  $3(x,y)  to  satisfy 
the  constraint  (7.2).   The  approximation  for  U(x,y)  was 

U(x,y)  =  Ci$i  +  C2$2  +  C3$3  (7.7) 

=  (x  -  x2)(y  -  y2)(Ci  +  xyC2  +  x2y2C3). 

The  values  of  Ci,  Cz9    and  C3  were  obtained  by  techniques 
described  in  Chapter  IV.   Computer  Output  4  shows  the  val- 
ues of  Ci,  C2,  C3  and  the  values  obtained  for  U(x,y)  at 
various  points.   In  this  method  it  was  found  necessary  to 
sub-divide  the  region  into  forty  equal  sub-intervals  for 
each  independent  variable  in  order  to  get  satisfactory 
accuracy. 

The  best  approximation  to  U(x,y)  was  obtained  by  Galer- 
kin's method  using  Eq.  (7«5),  where  $1  was  the  desired 
function.   There  was  no  error  and  the  method  was  able  to 
detect  that  this  was  the  case.   However,  when  Eq .  (7.6) 
was  used  the  maximum  error  was  eight  thousandths.  The  time 
required  to  solve  the  problem  by  this  method  was  0.57 
seconds . 

The  problem  was  solved  by  dynamic  programming  in  three 
seconds  with  accuracy  to  six  digits.   However,  over  half 
of  the  computer  time  was  spent  obtaining  inverses,  which 
could  have  been  fed  as  data  from  solving  Eq.  (6.1)  in  this 
particular  case. 
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The  Rayleigh-Ritz  method  required  just  under  ten  seconds 
and  had  a  maximum  error  of  two  thousandths. 

Because  of  the  time  required  and  the  accuracy  obtained 
by  Rayleigh-Ritz,  only  Galerkin's  method  and  dynamic  pro- 
gramming were  used  to  solve  Eq.  (7.3). 

Galerkin's  method  obtained  the  same  values  and  accuracy 
in  the  solution  of  Eq.  (7.3)  as  it  did  in  the  solution  to 
Eq.  (7.1)  and  took  the  same  time. 

Dynamic  programming  obtained  the  same  degree  of  accuracy 
in  solving  Eq.  (7.3)  as  it  did  in  solving  Eq.  (7.1)  and 
took  three  seconds.   The  results  are  shown  in  Computer 
Output  8  and  the  program  is  Computer  Program  7. 
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VIII.   CONCLUSION 

The  three  methods  considered  for  eigenvalue  problems 
yielded  satisfactory  values  for  the  eigenvalues  and  the 
eigenf unctions.   Generally,  the  relaxation  method  seemed 
to  be  most  satisfactory.   It  was  straight-forward  to  pro- 
gram, and  it  was  faster  then  Rayleigh-Ritz  and  dynamic 
programming.   It  converged  rapidly  even  if  step  functions 
were  used  on  several  different  tests  figures,  such  as  the 
L-shaped  and  triangular  regions . 

The  dynamic  programming  method  converged  in  the  same 
number  of  iterations  as  relaxation,  but  gave  poorer  esti- 
mates of  the  second  and  third  eigenvalues.   It  of  course 
was  much  more  difficult  to  program  and  required  more  com- 
puting time  due  to  the  needed  inverses. 

The  Rayleigh-Ritz  method  seemed  to  have  little  to  re- 
commend it  due  to  the  computer  time  required.   It  required 
much  finer  meshing  in  order  to  obtain  a  satisfactory  accu- 
racy.  The  only  advantage  it  had  was  that  no  iteration  was 
required . 

Of  the  three  methods  considered  in  the  solution  of 
Poisson's  and  the  biharmonic  equations   Galerkin's  method 
was  the  fastest  and  gave  accuracy  comparable  to  the  Ray- 
leigh-Ritz method.   It  was  also  the  simplest  of  the  three 
methods  used  to  program. 


49 


Dynamic  programming  gave  the  best  accuracy  generally, 
but  it  required  more  computer  time  than  Galerkin's  method 

Rayleigh-Ritz  had  the  same  difficulties  as  it  did  in 
the  eigenvalue  problem  and  was  considered  of  little  use. 

While  dynamic  programming  required  more  time  for  the 
solution  of  both  eigenvalue  problems  and  elliptic  partial 
differential  equations,  it  was  very  powerful.   It  only 
required  a  knowledge  of  the  function  U  on  the  boundary. 
It  can  be  extended  to  irregular  regions  £  1} .   It  obtained 
very  good  accuracy.   Much  of  the  time  was  spent  computing 
the  inverses.   If  the  same  points  were  used  in  solving 
several  different  problems,  these  inverses  could  be  cal- 
culated once  and  thereafter  entered  as  data,  reducing  the 
computer  time  greatly. 
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COMPUTER    OUTPUT    1 . ( REL AXATI ON ) 


EIGENFUNCTION=l 


TH 

=  8 

OMEGA= 

4,424858 

X 

Y 

U(X,Y) 

0 

• 

0 

• 

1 

0.09719044 

0 

• 

0 

• 

3 

0.2508801 

0 

• 

0 

• 

5 

0.3094832 

0 

• 

0 

• 

7 

0.2526796 

0 

• 

0 

• 

9 

0.09705645 

0 

• 

3 

0 

• 

1 

0.2508801 

0 

• 

3 

0 

• 

3 

0.6512997 

0 

• 

3 

0 

• 

5 

0.8061690 

0 

• 

3 

0 

• 

7 

0.6586339 

0 

• 

3 

0 

• 

9 

0. 2531579 

0 

• 

5 

0 

• 

1 

0.3094832 

0 

• 

5 

0 

• 

3 

0.8061690 

0 

• 

5 

0 

• 

5 

1.000000 

0 

• 

5 

0 

• 

7 

0.8169372 

0 

• 

5 

0 

• 

9 

0.3127754 

0 

• 

7 

0 

• 

1 

0.2526796 

0 

• 

7 

0 

• 

3 

0.6586339 

0 

• 

7 

0 

• 

5 

0.8169372 

0 

• 

7 

0 

• 

7 

0.6658412 

0 

• 

7 

0 

• 

9 

0.2551157 

0 

• 

9 

0 

• 

1 

0.09705645 

0 

• 

9 

0 

• 

3 

0.2531579 

0 

• 

9 

0 

• 

5 

0.3137754 

0 

• 

9 

0 

• 

7 

0.2551157 
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COMPUTER    OUTPUT    1 .( RELAX ATION ) 


EIGENFUNCTI0N=2 


TH= 

=  11 

OMEGA 

= 

6.927361 

X 

Y 

U(X,Y) 

0 

• 

0 

• 

1 

0.1853049 

0 

• 

0 

• 

3 

0.4886459 

0 

• 

0 

• 

5 

0.6086538 

0 

• 

0 

• 

7 

0.4958898 

0 

• 

0 

• 

9 

0.1895596 

0 

• 

3 

0 

• 

1 

0.3032387 

0 

• 

3 

0 

• 

3 

0.8008109 

0 

• 

3 

0 

• 

5 

0. 9981067 

0 

• 

3 

0 

• 

7 

0.8116993 

0 

• 

3 

0 

• 

9 

0.3099311 

0 

• 

5 

0 

• 

1 

0.001278793 

0 

• 

5 

0 

• 

3 

0.008393560 

0 

• 

5 

0 

• 

5 

0.01233941 

0 

• 

5 

0 

• 

7 

0.008128799 

0 

• 

5 

0 

• 

9 

0.001685064 

0 

• 

7 

0 

• 

1 

-0.3083511 

0 

• 

7 

0 

• 

3 

-0.8038315 

0 

• 

7 

0 

• 

5 

-0.9953710 

0 

• 

7 

0 

• 

7 

-0.8089499 

0 

• 

7 

0 

• 

9 

-0.3103257 

0 

• 

9 

0 

• 

1 

-0.1921543 

0 

• 

9 

0 

• 

3 

-0.50C6261 

0 

• 

5 

0 

• 

5 

-0.6194299 

0 

• 

9 

0 

• 

7 

-0.5025631 
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COMPUTER    OUTPUT    1  .  ( REL AXATION ) 


EI  GEN  FUNCTIONS 

PATH=11  01 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


2 

3 
3 
3 
3 
5 
5 
5 
5 
5 
7 
7 
7 
7 
7 
9 
9 
9 
9 


GA= 

6.927361 

Y 

U(X,Y) 

0 

• 

1 

0,1857997 

0 

• 

3 

0,3045443 

0 

• 

5 

0.002898388 

0 

• 

7 

-0. 3070407 

0 

• 

9 

-0.1916554 

0 

• 

1 

0.4894530 

0 

• 

3 

0.8029368 

0 

• 

5 

0.01102810 

0 

• 

7 

-0. 8017028 

0 

• 

9 

-0.4998162 

0 

• 

1 

0.6086553 

0 

• 

3 

0.9981196 

0 

• 

5 

0.01235584 

0 

• 

7 

-0.9953650 

0 

• 

9 

-0.6194320 

0 

• 

1 

0.4950783 

0 

• 

3 

0.8095817 

0 

• 

5 

0.005508117 

0 

• 

7 

-0. 8110784 

0 

• 

9 

-0.5033802 

0 

• 

1 

0. 1890559 

0 

• 

3 

0.3086166 

0 

• 

5 

5.9269C0«-05 

0 

• 

7 

-0.3116446 
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COMPUTER    OUTPUT    2 .( RELAXATION ) 


EIGENFUNCTION=l 
PATH=17  OMEGA= 


19.60767 


U(X,Y) 


0 

• 

1 

0 

• 

1 

0.1011795 

0 

• 

0 

• 

3 

0,2574610 

0 

• 

0 

• 

5 

0.3156426 

0 

• 

0 

• 

7 

C. 2600489 

0 

• 

0 

• 

9 

0.1008759 

0 

• 

3 

0 

• 

1 

0.2580318 

0 

• 

3 

0 

• 

3 

0.6597930 

0 

• 

3 

0 

• 

5 

0.8128256 

0 

• 

3 

0 

• 

7 

0.6714125 

0 

• 

3 

0 

• 

9 

0.2608653 

0 

• 

5 

0 

• 

1 

0.3149233 

0 

• 

5 

0 

• 

3 

0.8084638 

0 

• 

5 

0 

• 

5 

1.000000 

0 

• 

5 

0 

• 

7 

0.8271890 

0 

• 

5 

0 

• 

9 

0.3214785 

0 

• 

7 

0 

• 

1 

0.2583916 

0 

• 

7 

0 

• 

3 

0.6639176 

0 

• 

7 

0 

• 

5 

0.8214287 

0 

• 

7 

0 

• 

7 

0.6778069 

0 

• 

7 

0 

• 

9 

0.2627081 

0 

• 

<5 

0 

• 

1 

C. 1004035 

0 

• 

9 

0 

• 

3 

0.2579977 

0 

• 

9 

0 

• 

5 

0.3188558 
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COMPUTER    OUTPUT    2 .( RELAXATION  ) 


EI  GEN  FUNCTION- 2 

PATH=16  OMEGA= 


48.00555 


U(XTY) 


0 

• 

1 

0 

• 

1 

0.1909801 

0 

• 

0 

• 

3 

0.4904021 

0 

• 

0 

• 

5 

0. 6059132 

0 

• 

0 

• 

7 

0.4979762 

0 

• 

0 

• 

9 

0.1925782 

0 

• 

3 

0 

• 

1 

C. 3037586 

0 

• 

3 

0 

• 

3 

0.7884125 

0 

• 

3 

0 

• 

5 

0.9826701 

0 

• 

3 

0 

• 

7 

0.8096838 

0 

• 

3 

0 

• 

9 

0.3130895 

0 

• 

5 

0 

• 

1 

0.003429962 

0 

• 

5 

0 

• 

3 

0.01207077 

0 

• 

5 

0 

• 

5 

0.01897281 

0 

• 

5 

0 

• 

7 

0.01404283 

0 

• 

5 

0 

• 

9 

0.004551671 

0 

• 

7 

0 

• 

1 

-C. 3064085 

0 

• 

7 

0 

• 

3 

-0.7943833 

0 

• 

7 

0 

• 

5 

-0.9861545 

0 

• 

7 

0 

• 

7 

-0.8114642 

0 

• 

7 

0 

• 

9 

-0.3132687 

0 

• 

9 

0 

• 

1 

-0. 1937618 

0 

• 

9 

0 

• 

3 

-0.5022484 

0 

• 

9 

0 

• 

5 

-0.6227177 

55 


COMPUTER    OUTPUT    2 . ( REL AX ATION ) 


EIGENFUNCTI0M=3 
PATH=19  OMEGA= 


48.02386 


U(X,Y) 


0 

• 

1 

0 

• 

1 

0.1851881 

0 

• 

0 

• 

3 

0.30  20115 

0 

• 

0 

• 

5 

0.009329475 

0 

• 

0 

• 

7 

-0.3049613 

0 

• 

0 

• 

9 

-0.1935259 

0 

• 

3 

0 

• 

1 

0.4791670 

0 

• 

3 

0 

• 

3 

0.7890626 

0 

• 

3 

0 

• 

5 

0.02798434 

0 

• 

3 

0 

• 

7 

-0.7933792 

0 

• 

3 

0 

• 

9 

-0.5036968 

0 

• 

5 

0 

• 

1 

0.5902148 

0 

• 

5 

0 

• 

3 

0.9785894 

0 

• 

c 

0 

• 

5 

0.03599443 

0 

• 

5 

0 

• 

7 

-0.9824019 

0 

• 

5 

0 

• 

9 

-0.6229793 

0 

• 

7 

0 

• 

1 

0.4810330 

0 

• 

7 

0 

• 

3 

0.7979080 

0 

• 

7 

0 

• 

5 

0.02484581 

0 

• 

7 

0 

• 

7 

-0.80  38256 

0 

• 

7 

0 

• 

9 

-0.5073704 

0 

• 

9 

0 

• 

1 

0. 1851680 

0 

• 

9 

0 

• 

3 

0.3065956 

0 

• 

9 

0 

• 

5 

0.007537059 
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COMPUTER  OUTPUT  3.(GALERKIN) 


VALUES  OF  C(I)  ARE 


CI 

= 

1.000000 

C2 

= 

0 

C3 

= 

0 

U(X,Y) 


0 

• 

0 

• 

1 

0,008099992 

0 

• 

0 

• 

3 

0.01889999 

0 

• 

0 

• 

5 

0.02249999 

0 

• 

0 

• 

7 

0.01890001 

0 

• 

0 

• 

9 

0.008100022 

0 

• 

3 

0 

• 

1 

0.01889999 

0 

• 

3 

0 

• 

3 

0.04409999 

0 

• 

3 

0 

• 

5 

0.05249999 

0 

• 

3 

0 

• 

7 

0.04410003 

0 

• 

3 

0 

• 

9 

0.01890005 

0 

• 

5 

0 

• 

1 

0.02249999 

0 

• 

5 

0 

• 

3 

0.05249999 

0 

• 

5 

0 

• 

5 

0.06250000 

0 

• 

5 

0 

• 

7 

0.05250004 

0 

• 

5 

0 

• 

9 

0.02250007 

0 

• 

7 

0 

• 

1 

0.01890001 

0 

• 

7 

0 

• 

3 

0.04410003 

0 

• 

7 

0 

• 

5 

0.05250004 

0 

• 

7 

0 

• 

7 

0.04410006 

0 

• 

7 

0 

• 

9 

0. 01890007 
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COMPUTER  OUTPUT  ^.(RAYLEIGH) 


VALUES  OF    C(I)    ARE 

CI      =  1.120410 

C2      =  -0.5808935 

C3       =  0.3500372 


U(X,Y) 


0 

.25 

0 

.25 

0.03816108 

0 

.25 

0 

.50 

0.04937190 

0 

.25 

0 

.75 

0.03599290 

0 

.50 

0 

.25 

0.04937190 

0 

.50 

0 

.50 

0.06231650 

0 

.50 

0 

.75 

0.04461559 

0 

.75 

0 

.25 

0.03599290 

0 

.75 

0 

.50 

0.04461559 

0 

.75 

0 

.75 

0.03179573 
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COMPUTER  OUTPUT  5.(RAYLEIGH) 

EIGENFUNCTION=l  OMEGA=  4.447512 

THE    VALUES    OF    C( I )    ARE 

CI      =  29.89044 

C2       =         0.003808264 
C3       =  2.580980 

X  Y  U(X,Y) 


0 

.25 

0 

.25 

1.073552 

0 

.25 

0 

.50 

1.401157 

0 

.25 

0 

.75 

1.028184 

0 

.50 

0 

.25 

1.431359 

0 

.50 

0 

.50 

1.868153 

0 

.50 

0 

.75 

1.370868 

0 

.75 

0 

.25 

1.073485 

0 

.75 

0 

.50 

1.401070 

0 

.75 

0 

.75 

1.028116 
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COMPUTER    OUTPUT    5.(RAYLEIGH) 

EIGENFUNCTION=2  OMEGA=  7.230292 

THE    VALUES    OF    C(I)    ARE 

CI       =  9.664003 

C2      =  112.0784 

C3      =  -112.0294 

X  Y  U(X,Y) 


0 

.25 

0 

.25 

0.3401793 

0 

.25 

0 

.50 

1.766417 

0 

.25 

0 

.75 

2.309447 

0 

.50 

0 

.25 

-0.8598440 

0 

.50 

0 

.50 

0.6040002 

0 

.50 

0 

.75 

1.765845 

0 

.75 

0 

.25 

-1.629948 

0 

.75 

0 

.50 

-0.8604184 

0 

.75 

0 

.75 

0.3393202 

60 


COMPUTER  OUTPUT  5.(RAYLEIGH) 

EIGENFUNCTI0N=3  OMEGA=  7.247070 

THE   VALUES    OF   C( I )    ARE 

CI      =  9.656152 

C2      =  -112.4433 

C3      =  -111.6595 

X  Y  U(X.Y) 


0 

.25 

0 

.25 

-1.630179 

0 

.25 

0 

.50 

-0.8650630 

0 

.25 

0 

.75 

0.3325842 

0 

.50 

0 

.25 

-0.8558773 

0 

.50 

0 

.50 

0.6035087 

0 

.50 

0 

.75 

1.761141 

0 

.75 

0 

.25 

0.3463630 

0 

.75 

0 

.50 

1.770327 

0 

.75 

0 

.75 

2.309128 
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COMPUTER    OUTPUT    6. (DYNAMIC) 


U(X,Y) 


0 

• 

0 

• 

1 

0.008099854 

0 

• 

0 

• 

3 

0.01889967 

0 

• 

0 

• 

5 

0.02249958 

0 

• 

0 

• 

7 

0. 01889964 

0 

• 

0 

• 

9 

0,008099858 

0 

• 

3 

0 

• 

1 

0.01889966 

0 

• 

3 

0 

• 

3 

0.04409914 

0 

• 

3 

0 

• 

5 

0.05249893 

0 

• 

3 

0 

• 

7 

0.04409913 

0 

• 

3 

0 

« 

9 

0.01889965 

0 

• 

5 

0 

• 

1 

0.02249958 

0 

• 

5 

0 

• 

3 

0.05249896 

0 

• 

5 

0 

• 

5 

0.06249879 

0 

• 

5 

0 

• 

7 

0. C5249900 

0 

• 

5 

0 

• 

9 

0.02249960 

0 

• 

7 

0 

• 

1 

C.C1889965 

0 

• 

7 

0 

• 

3 

0.04409916 

0 

• 

7 

0 

• 

5 

0.05249896 

0 

• 

7 

0 

• 

7 

0.04409920 

0 

• 

7 

0 

• 

9 

0.01889972 

c 

• 

9 

0 

• 

1 

0. CC8099854 

0 

• 

9 

0 

• 

3 

0.01889966 

0 

• 

9 

0 

• 

5 

0.02249962 

0 

• 

9 

0 

• 

7 

0.01889972 

0 

• 

9 

0 

• 

9 

0.008099906 
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COMPUTER    OUTPUT    7. (DYNAMIC) 


EIGENFUNCT!ON=l 
PATH=    4  OMEGA= 


4.423473 


U(X,Y) 


0 

• 

0 

• 

1 

0.09554338 

0 

• 

0 

• 

3 

C. 2500933 

0 

• 

0 

• 

5 

C. 3091003 

0 

• 

0 

• 

7 

0.2500929 

0 

• 

0 

• 

9 

0.09554338 

0 

• 

3 

0 

• 

1 

0.2500930 

0 

• 

3 

0 

• 

3 

0.6546414 

0 

• 

3 

0 

• 

5 

0.8090985 

0 

• 

3 

0 

• 

7 

0.6546412 

0 

• 

3 

0 

• 

9 

0.2500929 

0 

• 

5 

0 

• 

1 

C. 3091002 

0 

• 

5 

0 

• 

3 

C. 8090987 

0 

• 

5 

0 

• 

5 

1.000000 

0 

• 

5 

0 

• 

7 

0.8090994 

0 

• 

5 

0 

• 

9 

C. 3091004 

0 

• 

7 

0 

• 

1 

0.2500930 

0 

• 

7 

0 

• 

3 

0.6546419 

0 

• 

7 

0 

• 

5 

0. 8090992 

0 

• 

7 

0 

• 

7 

0.6546427 

0 

• 

7 

0 

• 

9 

0.2500940 

0 

• 

9 

0 

• 

1 

0.09554315 

0 

• 

9 

0 

• 

3 

0.2500930 

0 

• 

9 

0 

• 

5 

0.3091C1C 
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COMPUTER    OUTPUT    7. (DYNAMIC) 


EIGENFUNCTI0N=2 
PATH=    6  OMEGA= 


6.927172 


U(X,Y) 


0 

• 

0 

• 

1 

0.1913275 

0 

• 

0 

• 

3 

0,5004815 

0 

• 

0 

• 

5 

0.6183150 

0 

• 

0 

• 

7 

0.5004799 

0 

• 

0 

• 

9 

0.1913269 

0 

• 

3 

0 

• 

1 

0.3092604 

0 

• 

3 

0 

• 

3 

0.8089828 

0 

• 

3 

0 

• 

5 

0.9994567 

0 

• 

3 

0 

• 

7 

0.8089805 

0 

• 

3 

0 

• 

9 

0.3092589 

0 

• 

5 

0 

• 

1 

1 

,959012'-06 

0 

• 

5 

0 

• 

3 

3 

,750642,-06 

0 

• 

5 

0 

• 

5 

2 

,407420,-06 

0 

• 

5 

0 

• 

7 

-4< 

,548319«-07 

0 

• 

5 

0 

• 

9 

-7 

.790408»-07 

0 

• 

7 

0 

• 

1 

-C. 3092566 

0 

• 

7 

0 

• 

3 

-0.8089775 

0 

• 

7 

0 

• 

5 

-0.9994535 

0 

• 

7 

0 

• 

7 

-0. 8089827 

0 

• 

7 

0 

• 

9 

-0.3092611 

0 

• 

9 

0 

• 

1 

-C. 1913257 

0 

• 

9 

0 

• 

3 

-0.5004784 

0 

• 

9 

0 

• 

5 

-0.6183146 
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COMPUTER  OUTPUT  7. (DYNAMIC) 


EIGENFUNCTION=3 
PATH  =  11  OMEGA* 


6.927173 


U(X,Y) 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


3 
3 
3 
3 
3 
5 
5 
5 
5 
5 
7 
7 
7 
7 
7 
9 
9 
9 


0 

• 

1 

0.1913257 

0 

• 

3 

C. 3092576 

0 

• 

5 

-1, 

►1992C4'-06 

0 

• 

7 

-0.3092595 

0 

• 

9 

-0.1913273 

0 

• 

1 

C. 5004785 

0 

• 

3 

0.8089781 

0 

• 

5 

-1, 

►801975«-06 

0 

• 

7 

-0.8089815 

0 

• 

9 

-0.5004808 

0 

• 

1 

0.6183135 

0 

• 

3 

0.9994553 

0 

• 

5 

1 

,738089'-06 

0 

• 

7 

-0.9994555 

0 

• 

9 

-0.6183147 

0 

• 

1 

0.5004803 

0 

• 

3 

C. 8089837 

0 

• 

5 

4 

.765801'-06 

0 

• 

7 

-0.8089792 

0 

• 

9 

-0.5004812 

0 

• 

1 

0.1913269 

0 

• 

3 

0. 3092602 

0 

• 

5 

2 

,622819»-06 
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COMPUTER    OUTPUT    8. (DYNAMIC) 


U(X,Y) 


0 

• 

0 

• 

1 

0. 008099731 

0 

• 

0 

• 

3 

0.01889934 

0 

• 

1 

0 

• 

5 

0.02249917 

0 

• 

0 

• 

7 

0.01889930 

0 

• 

0 

• 

9 

0.008099727 

0 

• 

3 

0 

• 

1 

0.01889932 

0 

• 

3 

0 

• 

3 

0.04409828 

0 

• 

3 

0 

• 

5 

0.05249788 

0 

• 

3 

0 

• 

7 

0.04409826 

0 

• 

3 

0 

• 

9 

0.01889930 

0 

• 

5 

0 

• 

1 

0.02249917 

0 

• 

5 

0 

• 

3 

0.05249792 

0 

• 

5 

0 

• 

5 

0.06249751 

0 

• 

5 

0 

• 

7 

0.05249795 

0 

• 

5 

0 

• 

9 

0.02249918 

0 

• 

7 

0 

• 

1 

0.01889931 

0 

• 

7 

0 

• 

3 

0.04409831 

0 

• 

7 

0 

• 

5 

0.05249793 

0 

• 

7 

0 

• 

7 

0.04409834 

0 

• 

7 

0 

• 

9 

0. C1889938 

0 

• 

9 

0 

• 

1 

0.008099716 

0 

• 

9 

0 

• 

3 

0.01889932 

0 

• 

9 

0 

• 

5 

0.02249920 

0 

• 

9 

0 

• 

7 

C. C1889938 

0 

• 

9 

0 

• 

9 

0.008099772 
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COMPUTER    OUTPUT    9. (DYNAMIC) 


EIGENFUNCTION=l 
PATH=    3  OMEGA= 


19.57639 


U(X,Y) 


0 

• 

0 

• 

1 

0,09549332 

0 

• 

0 

• 

3 

0.2500034 

0 

• 

0 

• 

5 

C. 3090197 

0 

• 

0 

• 

7 

0.2500030 

0 

• 

0 

• 

9 

0.09549332 

0 

• 

3 

0 

• 

1 

0.25C0032 

0 

• 

3 

0 

• 

3 

0.6545131 

0 

• 

3 

0 

• 

5 

0.8090189 

0 

• 

3 

0 

• 

7 

0.6545127 

0 

• 

3 

0 

• 

9 

0.2500032 

0 

• 

5 

0 

• 

1 

0.3090197 

0 

• 

5 

0 

• 

3 

0.8090194 

0 

• 

5 

0 

• 

5 

1.000000 

0 

• 

5 

0 

• 

7 

C. 8090203 

0 

• 

5 

0 

• 

9 

0.3090  200 

0 

• 

7 

0 

• 

1 

0.2500033 

0 

• 

7 

0 

• 

3 

0.6545139 

0 

• 

7 

0 

• 

5 

0.8090203 

0 

• 

7 

0 

• 

7 

0.6545147 

0 

• 

7 

0 

• 

9 

0.2500044 

0 

• 

9 

0 

• 

1 

0.09549326 

0 

• 

9 

0 

• 

3 

0.2500034 

0 

• 

9 

0 

• 

5 

0.3090  206 
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COMPUTER    OUTPUT    9. (DYNAMIC) 


EIGENFUNCTION=2 
PATH=    4  OMEGA= 


47.98485 


U(X,Y) 


0 

• 

0 

• 

1 

C. 1911072 

0 

• 

0 

• 

3 

0.5001713 

0 

• 

0 

• 

5 

0.6181207 

0 

• 

0 

• 

7 

C. 5001690 

0 

• 

0 

• 

9 

0.1911063 

0 

• 

3 

0 

• 

1 

C. 3091234 

0 

• 

3 

0 

• 

3 

0.8090394 

0 

• 

3 

0 

• 

5 

0.9998273 

0 

• 

3 

0 

• 

7 

C. 8090362 

0 

• 

3 

0 

• 

9 

0.3091207 

0 

• 

5 

0 

• 

1 

-9, 

,384C42»-07 

0 

• 

5 

0 

• 

3 

8, 

,968278«-06 

0 

• 

5 

0 

• 

5 

1 

,008165«-05 

0 

• 

0 

• 

7 

3, 

,054505«-06 

0 

• 

5 

0 

• 

9 

-4 

.620023,-06 

0 

• 

7 

0 

• 

1 

-C. 3091281 

0 

• 

7 

0 

• 

3 

-0.8090287 

0 

• 

7 

0 

• 

5 

-0.9998166 

0 

• 

7 

0 

• 

7 

-C. 8090367 

0 

• 

7 

0 

• 

9 

-0.3091323 

0 

• 

9 

.0 

• 

1 

-0. 1911175 

0 

• 

9 

0 

• 

3 

-0.5001757 

0 

• 

9 

0 

• 

5 

-0.6181255 
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COMPUTER    OUTPUT    9. (DYNAMIC) 

EIGENFUNCTION=3 

PATH=    7  OMEGA=  47.98486 


0 

• 

1 

0 

• 

1 

0 

• 

1 

0 

• 

1 

0 

• 

1 

0 

• 

3 

0 

• 

3 

0 

• 

3 

0 

• 

3 

0 

• 

3 

0 

• 

5 

0 

• 

5 

0 

• 

5 

0 

• 

5 

0 

• 

5 

0 

• 

7 

0 

• 

7 

0 

• 

7 

0 

• 

7 

0 

• 

7 

0 

• 

9 

0 

• 

9 

0 

• 

9 

U(X,Y) 


0 

• 

1 

0.1911126 

0 

• 

3 

C. 3091245 

0 

• 

5 

-9, 

,687110'-07 

0 

• 

7 

-0.3091258 

0 

• 

9 

-0.1911120 

0 

• 

1 

C. 5001729 

0 

• 

3 

0.8090308 

0 

• 

5 

-2, 

,217028«-06 

0 

• 

7 

-0.8090366 

0 

• 

9 

-0.5001741 

0 

• 

1 

0.6181221 

0 

• 

3 

0.9998204 

0 

• 

5 

1 

.343255'-06 

0 

• 

7 

-0. 9998225 

0 

• 

9 

-0.6181232 

0 

• 

1 

C. 5001751 

0 

• 

3 

0.8090382 

0 

• 

5 

5 

.371387'-06 

0 

• 

7 

-0.8090329 

0 

• 

9 

-0.5001737 

0 

• 

1 

0. 1911140 

0 

• 

3 

0.3091279 

0 

• 

5 

3 

.674680'-06 
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COMPUTER  OUTPUT  1  0.  ( GAL ERK IN ) 


VALUES  OP  C(I)  ARE 


CI 

= 

1,596382 

C2 

= 

1.596392 

C3 

= 

-2.598899 

U(X,Y) 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 

c 

0 
0 
0 
0 
0 
0 


3 
3 
3 

3 
3 
5 
5 
5 
5 
5 
7 
7 
7 
7 
7 


0 

. 

1 

0.002375632 

0 

• 

3 

0.01059511 

0 

. 

5 

0.01862749 

0 

. 

7 

0.02069907 

0 

. 

9 

0.01103618 

0 

. 

1 

0.01059507 

0 

. 

3 

0.03192532 

0 

. 

5 

0.04658196 

0 

• 

7 

0.04633235 

0 

. 

9 

0. 02294400 

0 

• 

1 

0.01862741 

0 

. 

3 

0.04658184 

0 

. 

5 

0.05916639 

c 

. 

7 

0.05281764 

0 

• 

9 

0.02397245 

0 

• 

1 

0.02069896 

0 

. 

3 

0.04633217 

0 

• 

5 

0.05281758 

0 

. 

7 

0.04240138 

0 

• 

9 

0.01732974 
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COMPUTER  PROGRAM  1. 
RELAXATION  METHOD  APPLIED  TO 
V2U  =  -X2U,  EQ.  (1.1) 

B'GIN 

COMMENT   FINDS   FIGfNVALUES    AND   EIGENVECTORS    BY 
RELAXATION    BASSD   DN    RAYL£IGH'S    FORMULA    nN    ANY 
TYp'::    OF    DOMAIN. 

FOP  MAL    PARAMETS  RS  : 

r,N\  HM'GA 

0M2  OMSGA    SQUARED 

OMC2        LAST    VALUE    OF    ONTGA    SQUAFED 

ZNDRM      MAX.    NORM 

Nl  ^AX.    NUMBER    OP    X    POSITIONS 

°l  MAX.    NUMBc*    0^     Y    POSITIONS 

HI  M^SH    OF    X    VADI  &BL  F 

H2  MH  SH    OF     Y    VA^IABLF 

N2  ARPAY    OF    LOW'R     POSITIONS    0-    Y 

N3  ARRAY    n  =    UDD"R     POSITIONS    OF    Y 

U»UlfZ    UNKNOWN   E I  GSNF UNCTIONS 

PH  POTENTIAL    =NjRGY 

K5  KINETIC    ENjRGY; 

RFAL    KE.PE.OM2tOM0  2.    R  A  ,  Z  S.OM.ZNQP  W; 
RrAL    H1,H2,C1 ,C2,C3,CATf DOG, CAT  1, DOG  1 ; 
INT  "Gr  R    M,  M,l.  ,Ml,Pl; 
R?AL    H,H3,H^; 
0M2 :=5C00.0; 
N  :  =  1 ; 

~RJE :  =0.002; 
L  :  =  1 ; 

CAT  1:  =  5000.0; 
R€ftD(NlfPl»HlfH2 ); 
H :=H1    ^2; 
H3:=HK  •■'?; 
H4:  =M2  —2; 
c*f  p  I  NJ 

RSAL    A50/ Y    X( 0: :N1); 
RIAL    ARRAY    Y(0::P1); 
I  NT   Gc  R    ARRAY    N2 »N3 ( D : : Nl  ) ; 
Rr"  AL    APR  AY    Z ,  LI .  U 1  (  0  :  :  N  1 . 0 : :  P 1 )  ; 
FnR    I :  =r.    u<N,T  TL    N1     ^3 
FOP    J:=0    UNTIL    PI    00 

B~GIN    UK  I.J  )  :=C.  •); 

Z(I ,J) :=C0; 

"END ; 
FOR     I:=0    UNTIL    mi    DO    X(I):=I^Hl; 
FCP    I:=0    UNTIL    PI    DO    Y(I):  =  I-VH2; 
FOR     I  : =0    UMTI  L     Ni    DO    R"  ADON ( N2 ( T ) ) ; 
PO  R    I  :  =0    U  NT  I L    N 1    00    R  •':•  A  DON  (  N  3  (  I  )  )  ; 
8  : 
IF    L=l    TH-N 

B^GIN    L  :  =  "»  ; 

-OR    T :=r    UNTI L    Nl     DO 

FOP    J:=C    UNTIL    PI    r;.n 

Z(I,J):  =  (X(I)-XU  ) ?    2 )    ( Y(J )-Y( J  P     2) ; 
r,n    yp    n ; 

IND; 

IF    L=2    THCN 
P.- GIN 

-OR     l:  =  0    "NTH    Nl    00 
FO^    J:=0    UNTIL    °l     DO 

BuGIM 

U( I.J ):=Z(I.J  ); 

Z(I  ,J)  :  =  (X(I )-X(I  )  -2)     (Y(J  )-Y(J  )**2)*(0.5-X(I)  ) 

END; 
CAT :=KH;L:=3;N:=1;0M2: =5000.0; 
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go  to  d; 

-fir:; 

IF    1=3    TM    M 
57GIN 

cuP     I:=0    "NTIL    Nl    Dn 
FOR    J:='     UNTIL    PI    DO 
BIG  I  "•' 

ill]  J  )  :  =  (X(  I  )-X(!  )     ^2)-  (Y(J  )-Y(  J  )  **2)  '  (0.5-Y(  J)  )  ; 

CAT  1:='k*  ;i_  :=4;Ms=  1;QM2:  =  5000.  0; 
GO    TO    ?; 


COMMENT    G;  t    s-    AND    *  E ; 
I:  P f : =k; :=0.0  ; 
-OF     I  :  =?    UNT7  l    Nl-1     DO 
FOR    J :  =  n! 2  (  I  )    UNTIL     N  3 (  I  )  -  1    DO 
Q-GIN 

Pi:=P~+( ((Z(I+3 ,J)-Z( I, J) )/Hl )**2 )*H 
+(( <Z( I,J+1)-Z< I. J) )/H2)  •     2P H; 
K-:=K:+(Z(  !,  J  )       2  )     Hi--  H2; 

"  MD ; 
CMT.  2:=0M?; 

PM2  :=Pr-/i<  f;  DM:  =  SORT(  T'?  )  ; 
Tc     (OM2>OM02)     OR     (N>40)     TH"N 

5 -GIN 

WR1 T~ (MN  IT    CONVERGING" ) ; 

WRI  T:  i  "N  =  "  )  ;  WP  I  TE  ON  (  N )  ; 

WRITr(,,0MP2    AMD    3^2    A  -  .  "  )  ;WR  I  T=DN  (  GM02,0W2)  ; 

on   to   R; 
=ND; 

IP     iRS  (  0M2-PM02  )<~R£    TH'-N 

IOCCMTROLC ?) : 

WRITS ("   "  );'•■'■"  IT?  ("   ");WPIT-("   »); 

WFIT'M"  COMPUTER  OUTPUT  1 .(  RELAXATION  ) 

wpt  t   (••  m);w«*it*!  ("  "  ) ; 

T  W  t  p  T  H  !    r\  C  T   7  '.    •  =  ]    • 

WRIT2("  £IGHNFUNCTION="tL-l) ;WRITE("    ")  ; 

INTCT  L0SIZ::=2; 

WRIT1' («  PATH  =  ,,,N," 

writ  Mn   M  )  ;WR  iTr(»  «  ) ; 

WR  I  Tc  (  "  X 

"  U  (  X  ,  Y )  «  )  ; 

WRIT;  (••    •'  )  ; 

=  0F    I:=l    ST™°    ?    UNTIL    M  -i    DO 
FOP.    J:  =  l    ST  •?    ?"    UNTIL    P1-]    On 
B'GIN    TMT  CI TL ^S I Z -: =9  ; 

WRIT?("  ",I    Div    Nlt"."»I    R5M    Nlt" 

J    DIV    P1,".»,J    R?M    Pl.Z(I.J)); 

WOTTT(M      ,1); 

*-  N1  D  * —  V  '  rl  * 

IF     (5iBS(OM2-OM02)<tRA)     AND    (L<4)     TH'N    GO    TO    6; 

IF    ABS(QM2-0M02)<lRA    TH5N    GO    TO    R; 

N:=N+1 ; 

COMMENT    G\:T    M~W    Z; 


OMSGA=",OM); 
Y  •• 


i 


>n 


I 


FOR     I : =1    UMT! L    Nl 
B-GIN 

FOP     J:=N2(I)+1    UNTIL    N3(I)-1    DO 
Z  ( I  ,  J )  :  =  (  H^  (  Z  ( I  + 1  ,  J  )  +  Z  (  I  -i  ,  J  )  )  +H3  •-  (  Z  ( I  t  J+  1  )  +  Z  (  I ,  J  - 1  )  )  ) 

/(2.*H3+2.    H4-0m2;>(H**2)  )  ; 
CND; 

L>2    TH:=N 
BrGIN 

CI :=D0G:=00G1 :=C.C ; 
PQR     T  :=i     "NTI  L    Nl-i     DO 
FOR    J:=M2(I)+1    UNTIL    N3(!)-l    DO 

sroiN 

C1:=CI+Z(I?J]  ■  m  ; 
OOG:=00G+Z( I, J )    M( I» J)    H; 
DOG! :=D0G1  +  Z( I  ,  J)    Ul  (  I, J  )    H; 
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C2:=D0G/CAT ;C3:=D0G1/CAT1 ; 

C0D    I  :=1    UNTI  L    Nl-1    DO 

FCR     J:=N2(T)+1    UNTIL     N3(I)-1    DO 

Z(7  ,J)  :  =  Z(ItJ)-Cl-C2>  U(  I, J  )-C3*UKI,JI  ; 

~ND; 
D: ZNOPM:=0.0 ; 
FOR    I:=C    UNTIL    Nl    DO 
FOR    J:=N2(!)     UNTIL    N3(U     DO 

BEGIN 

ZS:=ABS(Z( I, J  )  ); 

Ic    ZS>ZM3RM    TH5  N    ZNORV:=ZS; 

END; 

FHR    I:=0    UNTIL    Ml    DO 

FOP  J:=N2(I)  UNTIL  N3 ( I )  DO  Z ( I  , J )  : =Z ( I , J ) / ZNORM; 

GO  TO  A; 

END; 

R:.:  NO. 
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COMPUTER  PROGRAM  2. 
RELAXATION  METHOD  APPLIED  TO 
V^U  =  X2U,  EQ.  (1.3) 

B7GIN 

CPMMzNT    CINDS    cIG"NVALUES    AND    EIGHNVCTORS    BY 
RELAXATION    PASFD    ON    RAYLEIGH'S    FORMULA    FOR 
D      4(U)=W ■-■*?    U. 

«=OPMAL    PAR  A  MET  PS  : 

ON  OM=GA 

CM2  OM'^GA    SQUARED 

0MO2  L£ST    V*LUt    OF    OM-GA    SQUARED 

ZNORM         MAX.     NORM 

H  MZSH    SIZ5    OF    X    AND    Y 

UtUliZ      UNKNOWN   EIGENFUNCTIQNS 

PE  POTENTIAL    ENERGY 

K£  KINETIC    ENERGY; 

REAL    ARRAY    Z,U ,U1  (C  :  :  1 4 , 0  : :  1  4  )  ; 

REAL    KS . ?c  » OM  2»  OM0  2 » ERA  t  ZStOMtZNORM  f H .CAT.DOG .CAT1 • DOG1 ; 
REAL    C1tC2tC3; 
INTEGER    N,L; 
REAL    ARRAY    Y(0: :10); 
R-AL    ARRAY    X(0::10); 
LOGICAL    SWITCH; 
RC$L    CMn; 
OM: =7000.0; 
SWITCH  :=TRUE? 
CAT:=CAT1:=2GD0.0; 

H:=0.I;nM2:=2000.C;N:=l;  EP 4: =0.003 ;L:=1; 
OM 2:  =700000.0; 
t:  p  ft  j  =q    00  5* 

FOR* I  :=0   UNTIL    10    DO    X( I ) :=Y(I ) : *I*H; 

C0R    I:=0    UNTIL    14    D3 

ppp    J:=C    UNTIL    14    DO    Z ( I , J ) : =U ( I , J ) : =U 1 ( I , J ) : =0 . 0 ; 

BEGIN 

ppO'-'-nuR-   NnpM*;; 

BEGIN  ZNORM:-0.0; 
FHR  I: =3  UNTIL  11  ^n 
FOP    J: =3    UNTIL    11    03 

B'iGIN    ZS:=ABS(Z(I  ,  J  )  )  ; 

1=  ZS>ZNORM  THEN  ZNORM:=ZS;END; 
FOR  I:=0  UNTIL  14  DO 

FOR  J:=0  UNTIL  14  DO  Z ( I  , J  ) : =Z ( I , J ) / ZNORM; 
GO  TO  A; 

END  NORMA; 
PROCEDURE  PUNCT; 
BEGIN 
IF  L=l  THEN 

BEGIN  L:=2; 

C0R  l:=0  UNTIL  10  DO 

POP  J:=0  UNTIL  ID  DO 

Z(I+2,J+2):=(X(I )-X(I)**2)*<Y(J)-Y(J)**2); 

NORMA ;ENO; 
IF  L  =  2  TH5N 

BEGIN  L:=3;CAT:=<=;N:=l;0M2:=2000.0; 

CM: =7000.0 ; 

FOR  I:=0  UNTIL  14  DO  FOR  J : =0  UNTIL  14  DO  U ( I  ,  J  )  : =Z (  I , J  )  ; 

FOR  I:=0  UNTIL  14  DO  FOR  J : =0  UNTIL  14  DO  Z ( I  . J) : =0. 0; 

FOR  I:=0  UNTIL  ID  DO 

FOP  J:=0  UNTIL  10  DO 

Z(H-2tJ  +  2):  =  (X(  I>-X(I)**2)*(Y(J)-Y(J)**2)*(.5-X(I)  ); 

NOR  MA; END; 
IF  L=3  THEN 

BEGIN  L : =4 ; CAT l:=KE;N:=l;0M2:= 2000.0; 
OM:=7000.0; 

FOR  l:=0  UNTIL  14  DO  =0P.  J:=0  UNTIL  14  DO  Ul ( I . J) :=Z ( I . J) 
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cq»    I:=0    UNTIL    14    DO    FOR    J:=0    UNTIL    14    DO    Z ( I  ,  J ) : =0. 0  ; 

PHF     I :  =  0    UNTI L    10    DO 

FOR    J:=0    UNTIL     10    DO 

Z(I+2 ,J+2 ):  =  (X( I  )-X(I  )**2)*CY(J)-Y(J)**2)*<  .5-Y(J ) )  ; 

NOR  MA; "NO; 

END    FUNCT; 
Ic    SWITCH    THIN    BtGIN    SWITCH: =FALS 5 ; FUNCT  ;END; 
COMMENT    GCJT    P*    AND    K5  ; 
A: P1:=K~:=j.t  ; 
FOP    I :=2    UNTIL    12    DO 
FOR     J: =2    UNTIL    12    DO 

BEGIN    PJ:=^r+(Z( I +1 , J  )-4. *Z ( I , J ) +Z ( I- 1, J ) +Z ( I » J+ 1) 
+  Z(I  , J-l) )**2/<H*H); 

K/:=K-+(Z(  I, J  P-H)  **2;SND; 
0MO:=OM; 

QM2:  =  PS:/K?;0M:=SQRT(  0M2)  ; 
IF     (CM>OMO)    DR    (N>60)     THIN 

BIG  IN    WRIT:"  ("NOT    CTNV'RGING    "); 

WRIT=C    " ) ;WRITE( "NtZNORMfOMOfOM") ; 

WRIT£(N.ZNORM,OMO,OM) ;WRITE ("    " ) ; 

FOR    I :=1     UNTIL    11     DO 

WRITE(Z<3,I)»Z(5»I)  .ZC7.I)  .Z(9.I)  .Z(ll.I) ); 

GO    TO    Rs*ND; 
IF    ABS(C'MO-OM)<PRA    THEN 

BEGIN    IOCOMTROL (3)  ; 

WRITEC    '•);  WRITEC    ");WRITE(M    "); 

WRITEC  COMPUTER    OUTPUT    2.", 

"(RELAXATION) "  )  ; 

WRITEC"   ••) ;  WRITE (  "    ")  ; 

INTFISLDSIZE:=1; 

WRITEC  EIGENFUNCTI0N  =  ",L-1)  ;WR!TEC    "); 

INTFI5LDSI  ZE:=2; 

WRITEC  PATH=",N,"  OMFGA=" , CM ) ; 

WRITEC    ")  ; WRITEC    ")  ; 

WRITEC  X  Y  "» 

U(X.Y)")  ;  WRITEC    ");  WRITEC    "); 

PQR    I:=l    ST~P    2    UNTIL    9    DO 

FOR    J:=l     ST'.'P    2    UNTIL    9    DO 

BEGIN    INTFIf;LDSIZ::=2; 

WRITIC  ",I    DIV    10,".", I    REM    10, 

"  ",J    CIV    10,"."  ,  J    REM    10,Z(I+2, J+2) ); 

WRITS ("    " )  ; 
END; 

IF    L<4    THFN    FUNCT    ELSE    GO    TO    R; 
END; 
N:=N+l; 

COMMENT   GET   NEW   Z; 
FOR    I: =3    UNTIL    1 1    DO 
^OR    J:=3    UNTIL    11    DO 

Z(  I  ,J)  :  =  (  8.*Z<  1+1  ,  J)  +  8.*Z(I-1,  J)  +  8e*Z(I  ,J  +  1  ) 
+8.*Z(I, J-l)-Z(  T  +  2, J  )-Z(I,J-2)-2.*Z( 1+1, J  +  l  ) 
-2.*Z(I-1 ,J-1 )-2. - Z ( I -1 ,J+1  )-2.~Z(I+l,J-l  ) 
-Z(  1-2,  J  )-Z(I.J  +  2)  )  /(20.-QM2*H*H*H*H); 
FOR    I: =2    UNTIL    12    DO 

BEGIN    Z(0,U :=-Z(4,I) ; Z(14,I) :=-Z(10,I ); 
Z(1,I ):=-Z(3, I) ;Z( 13, I ) :=-Z( 11,1 ) ; 
Z(I  ,0) :=-Z(I,4) ;Z( I ,14):=-Z( I, IC  ); 
Z(I  ,1):=-Z(I,3);Z(I,13) :=-Z(I  ,11) ;END; 
IF    L>2    THEN 

BEGIN    CI :=D0G:=D0G1 :=0.0 ; 
FOR     I:  =3    UNTIL    11    DO 
FOR    J: =3    UNTIL    11    00 

BEGIN    Cl:=Cl+Z(I, J)*H*H; 
DOG:=DOG+Z( I, J )*U< I ,J)*H*H; 
DOGl:=OOGl+Z( I, J)*U1(I,J )«H*H;ENDJ 
C2:=D0G/CAT;C3:=D0G1/CAT1 ; 
FOR    I:=3    UNTIL     11    DO 
FOR    J:=3    UNTIL    11     00 

Z<I,J):=Z<IiJ)-Cl-C2*U<I,J)-C3*Ul<ItJ); 
FOR    I:  =  2    Uf'TIL    12    DO 
BIGIN 
Z(0,1  )  :=-Z(4,I  )  ;Z(  14,1)  :=-Z(10.I  );Z(1  .1  )  :=-Z(3.I  )*, 
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END 
IF    L<5 

END; 
P:FND. 


Z(13,I  ):— Zdli 
Z(I  ,1)  :=-Z(I,3) 

THEN    NORMA; 


T);Z(I,0):=-Z(I,4);Z(I,14) t=-Z(I ,10) 
; Z(  1,13): =-Z (1,2 1 );END; 
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COMPUTER  PROGRAM  3. 
GALERKIN'S  METHOD  APPLIED  TO 
V2U  =  2(x2+y2-x-y),  EQ .  (1.5) 


BEGIN 

CCMPENT    SOLVES    PARTIAL    DIFCERENTIAL    EQUATIONS    0- 
THE    FORM   L(U)=F    BY    GALERKIN'S    METHOD    ON    ANY    DOMAIN. 
FORMAL    PARAMETERS 


N 

M2 

M3 

N2 

N3 

F 

U 

C(I  ) 

Did) 

A(I  ) 
REAL 
REAL 
RfAL  H2 


NUMBER  OF  EQUATIONS 


POSITIONS 
POSITIONS 
POSITIONS 
POSITIONS 


MAX.  NUMBER  CF  X 

MAX.  NUMBER  OF  Y 

ARRAY  OF  LOWER  Y 

ARRAY  Oc  UPPER  Y 

KNOWN  FUNCTION 

UNKNOWN  FUNCTION 

ESTIMATING  FUNCTIONS 

^ESTIMATING  FUNCTIONS) 

COEF.  OF  ESTIMATING  FUNCTIONS 
H,E»OOG,CAT,OEV,AMUL,T, DOM; 
Hi; 
H3 


INTEGER  M2,M3; 
INTEGER  N,M,R ; 
E:=0.C00001;P:=0;DEV:=1.0; 
RSAD(N,M2»M3 ,H2,H3); 
Hl:=H2*H3; 


BEGIN 
REAL  ARRAY 


REAL 
REAL 

REAL 

REAL 
REAL 
INT  EG 


X(0: 
Y(0: 
3(1 


FCR 
FOR 
FOR 
FCR 
FOP 
FOR 


ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY  D, 
ER    AR^AY 


:M2); 
:M3); 
:Nfl: 

A  ,  C  ( 1 :  :  M ) 
F ,  U  (  0 :  :  M  2 


=0 
=0 
=0 
=0 

=  0 
=  0 


n+1  ) ; 

0  :  :  M  3  )  ; 
!1  (1  ::N  ,0  ::M2T0 
N2,N3( 0: :M2) 


M3  ) 


M2 
M3 
M2 
M2 
M2 
M3 


DO 
DO 
DO 
DO 
DO 
DO 


X(  I  ):  =  I*H2; 

Y(I):=I*H2; 
READONt N2(I ) ) 
READ0N(N3( I  )  ) 


UNTIL 
UNTI  L 
UNTIL 
UNTIL 
UNTIL 
UNTIL 
BEGIN 

F(I  ,J) :=0.0; 
U(I.J):=0.0; 
FOR    L:=l    UNTIL 

D(L,I ,J)  :=D1(L» I , J)  :=0.0; 
END; 
BEGIN    COMMENT    SOLVE 
PROCEDURE    WRITER ; 
BEGIN    WRITECSINGULAR"  ) 
END    WRITER; 
PROCEDURE    PERFORM; 


N 


HERE 


GO    TO    Q 


BEGIN    COMMENT 
WRI  TEC    ")  ;WRIT 


GET    U(  I 

(H       .1) 


") ;WRITE("    "  ) 


) 


WRI  TEC 

WRI  TEC 
WRI TEC  " 
WPITEC 
FCR    I :=1    UNTIL    N    DO 

BEGIN 

A(I ):=B(I.N+1) ; 

INTFIELDSIZE:=l; 

WRI  TEC 

END; 
FOR  J:=0  UNTIL  M3-1 
BEGIN 
FCR  I:=0  UNTIL  M2-1 

BEGIN  DOG:=O.C; 


J)  AND  A(  I  )  ; 
writs c  "); 

COMPUTER 


OUTPUT  3. (GALERKIN) " ) 


VALUES  OF  C(I)  ARE'MiWRITEC  "  ) 


C",I ,"=»,A(I) ); WRITE C  ") 


CO 
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CGP    L:  =  l    UNTIL    N   DO    DOG : =  DOG+ A( L )* D ( L  .  I . J ) ; 
U(I  ,J):=OQG; 
END; END; 

WRITEC"    ");WRITE("   ");WRITE<"    ") ; 
WRITEC  X  Y  ", 

"U(  x,Y)  ■•)  ; 
viritc  ("   ••)  ;WRIT5(M   "  ) ; 
FOR    !:=1    STEP    2    urjTI  L    M2-1    DO 
FOF     J:=l    ST=P    2    UNTIL     M3-1    DO 
R^C-IN    INTFIELDSIZE:=2; 
WRITEC  "»I    DIV    M2C'M    REM    M2, 

»  ",J    DIV    M3.".".J    REM    M3«U(ItJ)); 

WRITE (M    "  )  ; 
END  ; 
GC    TO    Q; 
END    PERcORM; 
PROCEDURE    SWITCH; 

BEGIN    COMMENT    SOLVES    FOR    A(I)    COEFF; 
INTEGER    Ml »Ml ; 
Ml:=N;Nl:=N+l; 

FOR    l:=l    UNTIL    Ml    DO     B(I •Ml+l ) :=C( I ) ; 
FOP    K:=l    UNTIL    Ml    DO 

BEGIN    IF    K=M1    TH=N    GO    TO    P; 
cOR    I:=K+1    UNTIL    Ml    DO 

BEGIN    IF    ABS(B(K,K)KABS(BUtKM    THEN 
BEGIN    R:=R+l; 
FOR    J:=l    UNTIL    Nl    DO 

BEGIN    T:=5(I,J);B(I»J):=B(K,J);B(K,J):=T;END; 
END;  END; 
P:IF    ABS(e(^tK))<F    THEN    WRITER    ELSE 

BEGIN    D£V:=DSV+ B(K,K); DCM:=B(K,K); 

FOP  J:=l  UNTIL  Nl  DO  B( K, J) : =B ( K , J ) /DOM; 

FOR  l:=l  UNTIL  Ml  CO 

BEGIN    AMUL:=B( I,K); 
IF    I  =  K    THEN       ELSE 
BEGIN 
FOR    J:=l     UNTIL    Nl     DO 

BEGIN    B( I ,JJ :=3(I . J ) -AMUL*B ( K  ,  J) ; 
END;END; END;£ND;END; 
P£RcOPM; 

END    SWITCH; 

CCMMENT  GUESSED  ^UNCTIONS  D(I)  HERE; 
FOP.  J:=C  UNTI  L  M2  DO 
FCP  K:=N2(J)  UNTIL  N3(J)  DO 
BEGIN 

F(  J  ,K)  :=2.*(X<J)**2  +  Y(KK*2-Y(K)-X(  J)  ); 
D(ltJ.K):=(X<J)-X(J)**2)*(Y(K)-Y(K  f**2)  ; 
C(2>J,K):  =  (X(J)**2-X(J)?"=3)SMY(K  )**.  2-Y  (  K  )+*  3  )  ; 
Cl(  I»J,K)  :=2.*<  X(  J)**2+Y(K)**2-X(  J)-Y(KH; 
D1(2.J.K):=(2.-6.*X(J) )*< Y(K )**2-Y(K )**3) 

+(2.-6.-Y(K) )t (X ( J)**2-X <J  )**3)  ; 
END; 
CCMKENT  SOLVE  cOR  INTEGRALS  F*D  AND  LD*D; 
FOR  I :=1  UNTI L  N  DO 
BEGIN  DOG: =0.0; 
FOR  J:=l  UNTIL  M2  DO 
FOP  K:=N2(J)+1  UNTIL  N3 ( J )  DO 
COG:=DOG+( F ( J ,K )* D( I , J ,K ) *H1 ) ; 
C(I  )  :=DOG; 
FOP  J: =1  UNTIL  N  DO 
BEGIN  CAT:=0.0; 
FOR  K:=l  UNTIL  M2  DO 
FOc  L:=N2(K)+1  UNTIL  N3(K)  DO 
CAT:=CAT+(D1( J,K,L)*D< I,K,L )*H1) ; 
B(I  ,J)  :=CAT; 
END; 
IF  I  =  N  THEN  SWITCH; 
END; 
END;  END; 
Q:END. 
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COMPUTER  PROGRAM  H. 

RAYLEIGH-RITZ  METHOD  APPLIED  TO 

V2U   =    2(x2+y2-x-y),    EQ.     (1.1) 


B-GIM 

CPMv~mt    SOLVTS    PARTIAL 

RAYL-IGH-PITZ    M«  THDO. 
FORV*L    PARAMETERS: 


DIcctR^NTIAL    EQUATIONS     BY 


H2  nx    S^CING 

H3  DY    SPACING 

N  NUMP-R     QP    1QU4TIONS 

F  KNOWN    FUNCTION 

^2  UPP~R    LIMIT    OF    X 

V3  UFP~R    LIMIT    CF    Y 

U  UNKNOWN   FUNCTION 

N3  UPP'R    BOUNDS    ON    Y 

N2  LOW-R    BOUNDS     DN    Y 

A(I)       CO'SF.    0^    ESTIMATING    FUNCTIONS 
"STTM/TTNn    FUNCTIONS; 
,DPG,CAT,ce V,AMUL,T, D0MTH1 ,H2  tH3; 


2  » 

C 


D(I  ) 

R?*L    H, 

int;7g-d  -.<h 

c  :  =0 .  00000 1 ;  R :  =0  ;  r\\r 
RFAD(N,  YI2,M3,H2,H3); 
HI:  =H2-^3; 

B7GTN 
RCAL    ARR£Y 

*RRAY 

ARRAY 

ARRAY 

ARRAY 

A  P  P  A  Y 


R?/a 

REAL 

RSAL 

R^AL 

RTAL 

INTcTG 

FOR    I 

FOP 

F0D 

cno 

FOP 
FOR 


X(C: 
Y(n: 

B(l 


M2); 

M3) ; 

:N,l: 

A  •  C  ( 1 : :  N  ) 

ctU (0:  :M2 
D,D1,D2(1 


R     ARRAY    N2»N3(C 


=0 
—  r\ 

=0 


=0 


M2 
M3 
M2 

M? 

M2 

M3 


Dn 
00 
DO 
DO 

DO 
DO 


UNTIL 
IJNTI  L 
UNTIL 

IJNJT  i_ 
UNTIL 
U  NT  I  L 

*:- GIN 

F(I  ,J)  :=U( I  ,J) 
=0C  L:=I  UNTIL 
D(L,I  ,J)  :=D1(L 

-MP  • 

FOR  "i: =1  UNTIL  N  DO 

FOR  J:=l  UNTIL  N+l  DO 

BJGTN 

COMMENT    SOIV.     H*R"; 

PROC-TURf    WRITR; 

Brr^TM    \-n  it":(  "S  TNGULAR"  ) 

INC    WRITER; 

ppCC-DURr    PERFORM; 

B"GIM 

g:;t  ij(i, 

")  ;WRIT£ 


=  1 


: N+l  )  ; 

Jo::M3) 

:  N ,  0  : : 
: :  M  2  ) ; 

X(  I  )  :=I 
Yd  )  :=I  <F3; 


M2  »C  :  :M3  ) 


0nN(M2(  I  )  )  ; 
ADON(  N3  (  I  )  )  ; 


=0.0; 

N    DO 

I, J) :=D2(L,I»J) :=0.0 


B(I TJ) 


•  c 


;o  to  o 


j) 


WRT   Tf  {  » 
W  R I T  -  ( •• 

WR7  T     (  "    tt)  ;  wo  IT;  ("    < 

WR  I  T" (  " 

C0R    I:=l    UNTIL    N    DO 

BEGIN 

A(! ):=B(I ,N+1) ; 

INT  FIELDS  I Z3:  =  l ; 

WR!  J.    (" 

*NC  ; 
FOR    I:=]    u  NT  T  |_    m  ?_ -I 
crp,    j:  =1     ! j>!t  t  i_    M3— f 

^VGI^    DOG: =0.0; 


AM  D    A  (  I  )  ; 
") ;  wri  t:  (•• 


•«) ; 

COMpUT" 

)  ;WRlTi("   •• ) ; 

VALU:S    OF    C(I) 


OUTPUT    A-.  (  RAYLf  I  OH  )  ")  ; 
AR=") ;WRITH(M    " ); 


!,"  =  ' 


(I)  );WRIT"("    ") 


DO 

CO 


:0P    L:=l    UNTIL    N    DO    DOG:  =  DOG+A(  L  )*D(  L  ,1  ,  J  ) 


79 


CIV 
100 


M2f".MT  1*100     Q~*''    M2, 


K  £  i 


*    M3,U( I t J) ) ; 


=C(I) ; 


U(I  ,J)  :=D^G; 

"ND; 
WRITE  ("    '•  );WRITr("    »  )  ;  WRI  T  i  ( '•    "); 

WR!  n  (  "  X 

"  U  (  X ,  Y )  "  ) ; 

writtC   ")  ;wpit^("   "  ) ; 
FOP    I:=0    ST~P    10    llNTTL    M2    00 
POP    J:=0    STEP     10    UNTIL    M3    00 
BEGIN    INTCI  ~LDSIZ::=2  ; 

wr  i  t:  (  »  »» ,  i 

writc("   »); 

~NC; 
on   to   0; 
END    PERFORM; 
PROCEDURE    SWITCH; 
B=GIN 

COMMENT    SHLV    S    FOP    ACT)     CC^F=; 
INT^GvR    NltMl; 
Ml: =N; Nl :=N+1 ; 

FOP    T:=l    UNTIL    Ml    DO    B(I,M1+1) 
FOP    K:=l    UNTIL    Ml    00 

BcGIN 

IP    k=mi    THEN   GC    TO    P; 

FOP     I:=K+1    UNTIL    Ml    00 
BEGIN 
IF    ABS(P( K,K) )<A3S(B( I ,«) )    TH~N 

BEGIN 

R  :  =P  +  i  • 

FOR    J:  =  l    UNTIL    Nl    00 

BEGIN    T:=8(I,J);3(I,J):=B(K,J) ;B ( K  ,  J ) : =T ;1ND; 
=ND;SMO; 

P:IF    A3S(3(K,K))<:     TH";N    WRITER    PLSE 
B'OIN 

D5V:=0HV*3(K, K) ;DOM:=B(K,K) ; 
=0R    J:  =  l    UNTIL    Nl    QP    B( K  .  J)  :=B ( K ,J) /DOM; 
-OR    I:=l    UNTIL    mi     DO 
B-GIN 

A'-'UL :=B(I,K); 
IP   i=k  th;n      *ls= 
b: "gin 

FOP    J:=!    UNTIL    Nl    DO 

r  :gin 

b(i  ,j)  :=b(i  ,j)-amul'p(k,j); 

;NC;  :ND;  =ND  ;H\'D;^NO; 
PERcORM; 
rND    SWITCH: 

COMMSNT    GU^SSiO    FUNCTIONS    0(1)    HEP"; 
Fnp     J:=0    UNT7L    M2    DO 
FOP    K:=N2(J)     UNTIL    N3(J)    00 

B?iGIN 

F(J,K)  :=2.   ■  <X( J)r    2+Y(K)*'2-X( J)-Y(K) ) ; 

D  <  1 .  J.K):=(X(J)-X(J)  **  2 )  •*(  Y  (  K )  -  Y(  K  )  **  2)  ; 

D(2,J,K):=X(J)    Y(K)-0(lvJfK); 

D(3,J,K):=X(J)    Y  (  * )    D ( 2  t  J , K ) ; 


=  (1.-2.     X(  J  )  )     (  Y(K  )-Y(K  )**  2)  ; 

=  (2.     X(  J)-3,  'X(J  )•-.   '■!)     (Y(K)v-2-Y(K)*c3); 

=  (  3  .     X  (  J  )  ■•  -.'  2  -4.     X  (  J  )  :  -  2  )      (  Y  (  K  )  >■  5  3  -  Y  (  K  )  ±*  4  )  ; 

=  (1.-2.    Y(K  )  )     ( X( J )-X( J  )'  ■  2)  ; 

=  (2,.     Y(K)-3.*Y(K)**.2  )'  (X(J)**2-X( J  )**  3  ); 

=  (3.     Y(K  )  ;  *  ?-4.   -Y(KH     3)     (  X(  J)    •    3-X(  .J)  -~  4)  ; 


OK  l.J.K  )  : 
DK2t  JfK): 
DIOtJtK)  : 
r2( 1,J,K  ): 

D2( ?fJ,K)  : 

D2(  3»J.k  ): 

"NO; 

CnMM"MT    SOLVz     ~^P,    INT    GR^LS    NOW; 
FOR     I:=l    UNTIL    N    DC 

B*GI  n 

DOG:=O.C ; 

FOP    J:=l    UNT7L    M2    DO 
=  0C    K:=M2(J)+1    UNTIL    N3(J)    DO 
OOr-  :=DOC— (  =(  J,K)*  D(I  , J,K)*H1)  ; 
C  ( I  ) :  =  DO  G ; 
rND; 
^OP    I  :  =1    UK'TT  L    N    CO 
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FPR    J:=l    UNTTL    N    DO 


7  M 


FO^    k  :  =  ]     'HTI  ' 
POP    L:=N2(K)+1 


M2    DO 
UNTIL 


N3(K)    DO 


B(!  fJ):  =  BtI,J)  +  (DKIf  KfL  )*Dl(J,K,L)+D2(IfK,L  )*D2( J,K,L)  )  * 


HI, 
7ND; 
SWITCH; 
=nd;*nd; 
Q:END. 
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COMPUTER  PROGRAM  5. 

RAYLEIGH-RITZ  METHOD  APPLIED  TO 

V2U  =  -A2U,  EQ.  (1.1) 

BEGIN 

COMMENT  SOLVES  EIGENVALUE  PROBLEMS  BY  RAYLEIGH-RITZ 
METHOO  USING  TRED2  AND  TGL2  FROM  NUMERISCHE 
MATHEMATIK  11, PP.  181-195  AND  PP.  293-306(1968)  BY  MARTIN 
REINSCH,BOWDLER»HILARY, AND  WILKINSON. 
FORMAL  PARAMETERS: 
N      NUMBER  OF  EIGENVALUES 

AND  THE  NUMBER  OF  EQUATIONS  USED 
HX     MESH  OF  X  VARIABLE 
HY     MESH  OF  Y  VARIABLE 
NPX    NUMBER  OF  X  POSITIONS 
NPY    NUMBER  OF  Y  POSITIONS 
U(I)   UNKNOWN  EIGENFUNCTIONS 
C(I)   ESTIMATING  FUNCTIONS  USED 
R(I)   COEF.  OF  ESTIMATING  FUNCTIONS 
D2(I)  EIGENVALUES  OF  MATRIX  A 
DKI)  EIGENVALUES  OF  MATRIX  B 
ZKI)  EIGENVECTORS  OF  MATRIX  B 
Z2(I)  EIGENVECTORS  OF  MATRIX  E=T'*A*T 
D(I)   USEO  IN  TRANSFORMATION  C=T«*D; 
INTEGER  NX,NY,N,NPX, NPY; 
REAL  HX,HY,EPS,DOG; 
LOGICAL  T; 
PROCEDUPE  TRED2( INTEGER  VALUE  N:REAL  ARRAY  D,E(*); 

REAL  ARRAY  A , Z ( *, * ) : RE AL  VALUE  PU ) : 
COMMENT  REDUCES  SYMMETRIC  MATRIX  TO  TRIDIG.  FORM  USING 
HOUSEHOLDER, S  REDUCTION. 
FORMAL  PARAMETERS: 

N   ORDER  OF  SYMMETRIC  MATRIX  A 
D   DIAGONAL  OF  RESULTS 
E   SUB-DIAGONAL  OF  RESULTS; 
BEGIN 

INTEGER  I,J,K,L; 
REAL  F,G,H,HH,TCL ; 
TOL:=PU/EPSILON; 
FOR  l:=l  UNTIL  N  DO 

FOR  J:=l  UNTIL  I  DO  Z (I  ,  J  )  : =A ( I , J ) : 
FOR  Il:=N  STEP  -1  UNTIL  2  DO 
BEGIN  I:  =  U; 

L:=I-2;F:=Z(I,I-1):G:=0.0; 
FOR  K:=l  UNTIL  L  DO  G: =G+Z( I , K ) *Z < I ,K ) ; 
H:=G+F*F; 
IF  G  <=  TOL  THEN 

BEGIN  E(I ) :=F  ;H:=0.0;GO  TO  SKIP;END: 
L:=L+1 : 

G:=E(I):=IF  F>=  0.0  THEN  -SQRT(H)  ELSE  SQRT(H): 
H:=H-F*G;Z( 1,1-1)  :  =  F-G : F: =0 .0 ; 
FOR  J:=l  UNTIL  L  DO 
BEGIN 

Z( J, I ) :=Z(I ,J)/H:G:=0.0; 

FOR  K:=l  UNTIL  J  DO  G: =G  +  Z ( J , K) *Z ( I , K  )  ; 
FOR  K:=J+1  UNTIL  L  DO  G: =G+Z ( K, J ) *Z ( I , K ) ; 
E( J) :=G/H:F:=F+G*Z(J,I) ; 
END; 

HH:=F/(H+H) : 

FOR  J:=l  UNTIL  L  DO 

BEGIN 

F:=Z(I,J) ;G:=E(J) :=E( J)-HH*F; 

FOR  K:  =  l  UNTIL  J  00 

Z( J,K)  :=Z( J  ,K)-F*E(K)-G#Z( I,K) : 

END; 
SKIP:D(I ) :=H; 
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G:=G+Z( I,K)*Z(K, J ) : 
Z(K,J):=Z(K, J)-G*Z(K, I) ; 


TRANSFORMATION; 


=  E(  I) 


END; 
D(l)  ;=E(l) :=0.0; 
FOR  l;=l  UNTIL  N  DO 
BEGIN  L:=I-l; 
IF  D(  I  )  -.=  0.0  THEN 
FOR  J:=l  UNTI L  L  DO 
BEGIN  G;=0.0; 
FOR  K;=l  UNTIL  L  DO 
FOR  K:=l  UNTIL  L  DO 
END; 
D(I ) ;  =  Z(I ,1) ;Z(I t I  ):  =  1.0; 

FOR  J;=l  UNTIL  L  DO  Z ( I  ,  J ) : =Z ( J  ,  I  )  : =0.0 ; 
END; 
END  TRED2: 
PROCEDURE  TQL2(INTEGER  VALUE  N;REAL  ARRAY  D,E(*) 

REAL  ARRAY  Z( * ,*)  J  LOG ICAL  RESULT  FAIL); 
COMMENT  FINDS  EIGENVALUES  AND  EIGENVECTORS. 
FCRMAL  PARAMETERS; 

N   ORDER  TRIDIAGONAL  MATRIX 
D   DIAGONAL  ELEMENTS 
E   SUB-DIAGONAL  ELEMENTS 
Z   MATRIX  OF  HOUSEHOLDER 
BEGIN 

INTEGER  I,J,K,L,M: 
REAL  B,CtF,G,H,P, R,S; 
FAIL:=FALSE; 

FOR  I :=2  UNTIL  N  DO  E( 1-1) 
E(N)  ;=B;  =  F;=0.0; 
FOR  Ll;=l  UNTIL  N  DO 
BEGIN  L:=Ll;J:=0; 

H:=EPSILON*(ABS(D(L) ) +ABS ( E ( L ) ) ) ; 
IF  B<H  THEN  B :=H; 

COMMENT  LOOK  FOR  SMALL  SUB-DIAGONAL  ELEMENT; 
FOR  Ml  ;=L  UNTIL  N  DO 
BEGIN  M:=M1 ; 

IF  ABS(E(M))<=B  THEN  GO  TO  CONT;END; 
CONT;IF  M=L  THEN  GO  TO  ROOT; 
NEXTIT;IF  J=30  THEN 

begin  fail;=true;go  to  fin;END; 
J;=J+1 ; 

P:=(D(L+1)-D( LI )/(2.*E(L) ) : 
R:=SQRT(P**2+1.0)  ; 

H;=D(L)-E(L)/( IF  P<0.0  THEN  P-R  ELSE  P+R); 
FOR  I:=L  UNTIL  N  DO  D ( I ) : =D( I ) -H; 
F ;  =F+H  * 

COMMENT  QL  TRANSFORMATION; 
P;  =  D(M) :C:  =  l.o;S;=O.C  ; 
FOR  I:=M-1  STEP  -1  UNTIL  L  DO 

BEGIN 

G;=C*E( I) ;h:=c*p; 

IF  ABS(P)  >  =  ABS(E( I ) )  THEN 
BEGIN 

C;=E(  I ) /P:R:  =  SQRT(C**2+1.0) ; 
E(  1  +  1 )  :=S*P*R;S:=C/R:C:=1.0/R: 
END  ELSE 
BEGIN 

C:  =  P/E(  I  )  ;R:  =  SQRT(C**2+1.0) ; 
E(I+1):=S*E(I )*R;S:=1.0/R;C:=C/R; 
END; 
P:=C*D(I)-S*G; 
D( 1+1 ) :=H+S*(C*G+S*D( I) ) ; 
COMMENT  FORM  VECTOR; 
FOR  K;=l  UNTIL  N  DO 
BEGIN 

H:=Z(K, 1+1) :Z(K,I+1):=S*Z(K,I)+C*H; 
Z(K,I ) :=C*Z(K, I )-S*H: 
END; END; 
E(L) ;=S*P;D(L ) ;=C*P; 

IF  ABS(E(D)  >  B  THEN  GO  TO  NEXTIT; 
ROOT;D(L) :=D(L)+F; 
END: 
COMMENT  ORDER  EIGENVALUES  AND  EIGENVECTORS: 
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FOR  l:=l  UNTIL  N 
BEGIN  K:=I:p: 
FOR  J: =1+1  UN 
IF  D(J)  <  P  T 
BEGIN  K:= 
IF  K  -.=  I  THE 
BEGIN 
D(K) :=D(I 
FOP  J:=l 
BEGIN  P 
END:END; 


DO 

=D(I ) : 

TIL  N  DO 

HEN 

J;P:=D(J) ;END; 

N 

)  :D( I ) :  =  P; 
UNTIL  N  DO 
:=Z( J, I)  :Z(J,I):=Z(J,K);Z(JfK) 


=  P;END: 


FIN 

END  TQL2; 
START:READ(N,NPX, 
BEGIN 

COMMENT  PROGRAM  STAR 
REAL  ARRAY  Dl ,D2,E1, 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 


NPY,NX,NY,HX,HY,EPS); 


REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 


TS  HERE; 
E2(l: :N) : 
1: :N)  : 

:NPX,0: :NPY) ; 
: :N,1  ::N)  ; 


D,R(1  : :N, 

U( 1 ::N,0: 

E,A,B,P(  1 

X(u: :NPX) 

Y(0:  :NPY) 

C,C1,C2(1 

Z1,Z2(1 
PROCEDURE  FCN; 
COMVENT  COMPUTES  FUNCTIONS  AND  MATRICES  A  AND  B: 
BEGIN 
X(O)  :=Y(0) :=C.O; 


: :N,0: :NPX,0::NPY) ; 
N, 1 : : N  )  ; 


FOR 
FOR 
FOR 
FOR 
FOR 
FOR 


=  1 
=  1 
=  1 
=  1 
=  0 
=0 


UNTIL 
UNTIL 
UNTIL 
UNTIL 
UNTIL 
UNTIL 


NPX  D 
NPY  D 
N  DO 
N  DO 
NPX  D 
NPY  D 


Cl( 


BEGIN 

C(1,I,J) :=(X( I )-X 
C(2,  I  ,J) :  =  (C.5-X( 
C(3tI ,J) :=(0.5-Y( 
J) :=(l.-2. 
=(Y( J)-Y(J 
J) :=(l.-2. 
J) :=(l.-2. 
J) :=(0.5-X 
J) :=(X( I )- 


It  I 
f  J) 


3,1 
1,1 
2,1 
3,1 


UNTIL 
UNTIL 


DO 
DO 


Cl( 
2,1 
Cl( 
C2( 
C2( 
C2( 
END 
FOR  K:=l 
FOR  L:=l 
BEGIN 
FOR  I:=l  UNTIL  NP 
FOR  J:=l  UNTIL  NP 
BEGIN 

A(K,L) :=A(K,L 
+  C2(K,I  ,  J)* 
B(K,L)  :=B(K,L 
END: 
END; 
END  FCN: 

PROCEDURE  WRITER: 
COMMENT  WRITES  ALL 
BEGIN 

UNTIL  N  DO 
UNTIL  I  DO 
D0G:=Z1( I, J 
UNTIL  N  DO 
UNTIL  N  DO 
UNTIL  N  DO 
UNTIL  N  DO 
DOG:=0.0: 
1  UNTIL  N 
=  DOG; 


0  X( I ) :=I*HX: 
0  Y( I ) :=I*HY; 

A( I, J) :=B( I , J) :=0.0; 

0 
0 

( I )**2)*(Y( J)-Y( J)**2) : 

I ) )*C<1,  I,  J)  : 

J)  )*C(1, 1  ,J)  ; 

*X(I  ) )*( Y( J)-Y( J)**2) ; 

)**2)*(0.5-3.*X( I )+3.*X( I )**2) : 

*X(I )  )*(0.5-Y( J)  )*(Y( J)-Y( J  )**2) : 

*Y< J ) )*( x< I )-x< I )**2) : 

( I ) )*(X( I )-X( I )**2)*( l.-2.*Y(J)); 

X( I )**2)*(U.5-3.*Y( J)+3.*Y( J  )**2) 


FOR  I:=l 
FOR  J:=l 
BEGIN 
FOR  L:=l 
FOR  l:=l 
FOR  L:=l 
FOR  l:=l 
BEGIN 
FOR  J: 
R(I ,L) 
END; 
FOR  l:=l 


X    DO 
Y    DO 

)+Cl (K,I , J)*C1 (L»I , J)*HX*HY 

C2(L, I, J)-HX*HY; 

)+C(K,I , J)*C(L,  I, J)*HX*HY: 


data; 

);Z1(I,J):=Z1(J,I):Z1(J,I) :=COG; END; 
D(L,I ) :=Z2( I,L) ; 

DO    D0G:  =  D0G+Z1U  ,J)*D(L,J)  : 


UNTIL    N    DO 
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BEGIN 


J:=0  UNTI 
L:=0  UNTI 
BEGIN  U( I 
FOR  K:=l 

C(K, J,L 
;END; 
FOR  K:=l  UNTIL  N 
BEGIN  IOCONTR 


FOR 
FOR 


END 


WRITEC 
WRITEC" 

WRITEC 
WRITEC 
WRITEC 
WRITEC 
WRI  TEC 
FOR 


")  ;wr 

M)  ;wr 

" )  :  wr 

" )  ;WR 
L:=l  UNTI 
BEGIN  INT 
WRITEC 
WRITEC  " 
WRITEC 
WRI  TEC  ")  ;WR 
FOR  J:=0  STEP 
FOR  L:=0  STEP 
BEGIN  INT 
WRITE(  " 

L  DIV 
WRITEC  " 
END;END: 
GO  TO  Q; 
END  WRITER; 
FCN: 

TRED2CN,D1,E1  ,B, 
TQL2CN,D1,E1,Z1, 
IF  T  THEN  BEGIN 
BEGIN 
FOR  I 


OMEGA=",SQRT(D2(K) ) ) 


U  ( X ,  Y )  " ) ; 


=1  UNTI 
=1  UNTI 
=1     UNTI 


FOR  I 
FOR  J 

end; 

FOR  l;=l  UNTIL  N 
FOR  K:=l  UNTIL  N 

BEGIN  P( I ,K) ; 

FOR  J:=l  UNTI 

end; 
for  i :=1  until  n 

FOR  J:=l  UNTIL  I 
BEGIN  DOG:=Zl 

FOR  l:=l  UNTIL  N 

FOR  K;=l  UNTIL  N 
BEGIN  E( I »K) : 
FOR  J:  =  l  UNTI 
END; 

TRED2(N,D2,E2,E, 

TQL2(N,D2,E2,Z2, 

IF  T  THEN  BEGIN 

WRITER; 

END; 

Q:G0  TO  START; 

END, 


L  NPX  DO 

L  NPY  DO 

,JtL) :=0.0; 

UNTIL  N  DO  U(I,J,L)  :=UCI,J,U+RC  I,K)* 

)  : 

DO 
0L(3)  ; 

ITE(  »«  ")  ;  WRI  TEC  "  ) ; 

COMPUTER  OUTPUT  5. ( RAYLE  IGH ) " ) ; 

ITEC  ")  ;WRITEC  "); 

EIGENFUNCTION=",K," 

ITE(  "  ")  tWRITEC  ••)  ; 

THE  VALUES  OF  C(I )"); 

I  TEC"  ") ; 

L  N  DO 
FIELDSIZE;=l; 

C",L ,"=",RCK,L) ) ; 
);WRITEC  ")  ;END;WRITE("  "); 

X  Y 

ITEC  ")  ; 
nx  until  npx  do 
ny  until  npy  do 
fieldsize;=2; 

",J  DIV  NPX,".", J  REM  NPX," 
NPY,".",L  REM  NPY,U(K,J,L) )! 

) ; 


zi,EPS) ; 

T)  : 

writecfail=1")  ;go  to  q;end  else 

L  N  DO  Did)  ;=SQRT(D1  (  I  )  ); 

L  N  DO 

L  N  DO  Zl(  J, I  )  :=Z1(  J,  II/OHI)  ; 

DO 

DO 

L  N  DO  P(I,K):=P(I,K)+A(I,J)*Z1(J,K): 

DO 

DO 
(I,J);Z1(I,J):=Z1(J,I);Z1(J,I) := DOG; END; 

DO 

DO 
=  0.0; 
L  N  DO  EU,K)  :=E(I,K)+Z1(I,J)*P(J,K): 

Z2,eps)  ; 

T)  : 

write("FAIL=2") ;G0  to  q?end; 
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COMPUTER   PROGRAM    6. 
DYNAMIC    PROGRAMMING    APPLIED   TO 


(1) 


V2U 


=    2(x2+y2-x-y) 


EQ.     (1.5) 


(2)    VU    =    A2U   EQ.     (1.3) 

BEGIN 

CCMVSNT    SOLVES    PARTIAL    DIFFERENTIAL 
F I  GEN  VALUE    PROBLEMS    ON    A    RECTANGLE 
PROGRAMMING    AND    STODOLA'S    METHOD. 
FORMAL       PARAMETERS: 
OM         OMEGA 
0M2      OMEGA    SQUARED 
N  NUMBER    OF    X    POSITIONS 

M  NUMBER    OF    Y    POSITIONS 

Kl         COEFF    IF    NEEDED 
VI         ORDER    OF    THE    OPERATOR 

V1=0      FIRST    ORDER    EIGENVALUE    PRB 
Vl=2       SECOND    ORDER    EIGENVALUE    PRB 
Tl  TYPE    OF    PROGRAM    RUN 

T1=0      FIRST    ORDER 
Tl=2      SECOND    ORDER 
Tl=50    EIGENVALUE 
P(I)    SECOND    NORMAL    DERIVATIVE    ON    BDRY 
El  FUNCTION    TO    SATISFY 

U  UNKNOWN    FUNCTION 

CI         NUMBER    OF    EIGENVALUES; 
REAL    KEjPE; 

REAL  e»ERA,OM,OM2,OMD2,H,Kl ; 
REAL  CAT1, CAT2.CL.C2tC3,DOGl,DOG2; 
REAL  ZNORMtZS; 
INTEGER  N,M,N1,S,T1, VI ,M1; 
INTEGER  01,02; 
LOGICAL  TOGGLE; 
REAL  ZNORMl; 

X:RSAD(N,M,H,Kl,Tl,QlfVl); 
TOGGLL :=TRUE; 
0  2 :  =  1 ; 

CAT1:=CAT2:=50C0.0; 
Ml:=2*M-2; 
0M2: =1000.0; 
Nl:=0; 
ERA:=0.03; 
ERA:=0.002; 
E:=0. 000001; 
BEGIN 

REAL  ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
A  p  c<  B  v 

ARRAY  I 1,0(  l:  :M-1.1: :M-1) 

ARRAY  U,8tR(0: :N,0: :M) 

ARRAY  A( l: :N,l: :M,1 : :M) ; 

ARRAY  D( l: :N.O: :M.O: :M1) 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

=0 

=0 


fOUATIONS  AND 
BY  DYNAMIC 


REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

R^AL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

FOR 

FOR 

FOR 

FCR 

FOR 

FOR 

IF 


N); 
M); 
N ,  0 : 

r:  :N 

N,0: 


Xl(  0: 
Y  1  ( 0  : 
U5(0: 
U1.U2( 
U3  ( 0 : 

CI  ,B1,P1( 0 
Z(0  ::N+4,0: 
El (0: :N,0:: 
I1,0(  l:  :M-1 
U  ,  B  T  R  ( 0 :  : 
A  ( 1 : :  N  f  1 : 
D( l: :N.O: 
CAT(0: :M) 
C  (  0  :  :  N  ,  0 : 
P 1  *  P  3  ( 0  :  : 
P2  .  PMO:: 


M) 


0: 

M) ; 

:  N  ,  0 : 
:M+4) 
f); 


M) 


:M); 

M)  ; 

N)  ; 
I:=0  UNTIL  M  DO  READOM(U(0»I ) ) ; 
I:=0  UNTIL  N  DO  READON(U(I,0) ) ; 
I:=0  UNTTL  M  DO  READON(U( N,I ) ) ; 
I:=0  UNTIL  N  DO  RE ADOM( U( I ,M ) ) ; 
1:=Q    UNTIL  N  DO  X1(I):=I*H; 
J:=0  UNTIL  M  DO  Y1(J):=J*H; 
(T1>0)  AND  (Vl=2)  THEN 
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=0    UNTIL    M    DO 

=  0.0; 

=0    UNTIL    M    DO 

=-PK  I)  ; 

=0    UNTIL    N    DO 

—  **)   r<  • 

=0*  UNTIL    M   DO 
=-P2(I) ; 


PEG  IN 
C0R  I 
DKI  ) 
FOR  I 
P3(  I  ) 
FOP  I 
P2(I  ) 
FOR  I 
P4(I  ) 

eNn ; 
FPP~I:=0    UNTIL    N    DO 
FOR    J:=0    UNTIL    M    DO    Ul ( I , J ) : =U2( I, J ) : =U3( I , J ) : =U5( I , J ) : =0. 0; 
S:=N; 
BEGIN 

PROCEDURE    WRITER; 
BEGIN    WPITE(  "SINGULAR")  ; 
GO   TO   V»; 
END    WRITER; 
PROCEDURE    WPITl; 
BEGIN 

IOCCNTROLO)  ; 

WRIT£(«    ");  WRITE  ("    ");  WRITE  ('•    ")  ; 

WRITEC  COMPUTER    OUTPUT    9.  (  DYNAMIC  )") ; 

WRITEC    ");  WRITEC    ");  WRITEC    "); 
INTFIELDSIZE:=1; 

WRITEC  EIGENFUNCTIPN  =  ",Q2-1)  ;WRITEC    " ) ; 

INTFISLDSI  ZE:=2; 

WRITEC  PATH  =  "»NK" 

WRITEC    '•)  ;WRITE("   "  )  ; 


OMEGA=".OM) ; 
Y  », 


DO 

DO 

",I     DIV    N«M."*I    REM    N." 


STARP  ; 


WRI TE (  " 

U(XtY)"); 

WRITEC    ");WRITE("    "  )  ; 
FOR    I :=1    STEP    2    UNTI L    N-l 
FOR    J:=1    STEP    2    UNTIL    M-l 
BEGIN    INTFIELDSIZ£:=2 ; 
WRI  TEC 

J    DIV    M,».»,J    REM    M,U(  It  J)  )  ; 
WRITEC    ") ; 
END; 
IF    Q2  =  01+l    THEN    GO    TO    W    ELS; 
END    WPITl; 
PROCEDURE    STARP; 
BEGIN 

IF    Tl=50    THEN 
BEGIN 
IF    02=1    THEN 

BEGIN    Q2:=2; 
FOR    I:=l    UNTIL    N-l    DO 
^OR    J:=l     UNTIL    M-l    DO 

U(  I  ,J)  :=X1(  I  )*Y1(  J)-(  l.-XKI)  »*  (l.-Yl(J)  ); 
GO   TO    Z5; 
END; 
IF    Q2=2    THEN 
BEGIN 

FOR    I :=1    UNTI L    N-l    DO 
FOR    J:=l    UNTIL     M-l    DO 
BEGIN 

U2(I tJ) :=U(Ii J); 
U5(Ii J ):=0.0; 
U(I  ,J):=X1  ( I  )•"•(  l.-XKI  )  )*(.5-XK  I)  )*Y1(J  )^(  l.-YKJ)  )  ; 

END; 
CATl:  =K E ; Q2 :  =  3 ;Nl:=l;0M2: =5000.0; 
GO    TO    Z5; 
END; 
IF    Q2=3    THEN 
BEGIN 

FOR    I :=1    UNTI L    N-l    DO 
FOR    J:=l    UNTIL    M-l    DO 
BEGIN 

U3(  K J) :=U(I , J)  ; 
U5(I » J):=0.0; 
U(I  ,J)  :=YK  J)*  (1.-Y1  (  J)  )*(.5-Yl  (  J)  )*X1(  I  )*<  l.-XKI  )  ); 
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END; 
CAT2:=K5;Q2:=4;0M2:=5000.0; 

GO   TO    Z5; 

END;  END; 
Z5: 

IF     Q2>2    THEN    GO    TO    Z6; 
COMMENT    DUT    FUNCTION    TO    SAT.    HERE; 
IF    Tl-«=50    THEN 

BEGIN 

FOR     I  :=0    UNTIL    N    DO 
FOR    J:=0    UNTIL    M    DO 
BEGIN 

E1(I,J)  :  =  2.*(  XHI  )**2+YK  J)**2-X1(I  )-YH  J)  )  ; 
C(It J):=2.0*H*H*£1(I,J)  ; 

END; END; 
FOR    l:=l    UNTIL    M-l    DO 
FOR    J:=l    UNTIL    M-l    DO 

BEGIN 

IF    I  =  J    THEN    I1(I.J):=1.0    ELSE    1 1 ( I  .  J ) :=0.0; 

Q  ( I  j  J  )  :  =0  . 0  ; 

IF    I=J    THEN    Q(I,J):=Kl*3.0; 

IF    ABS(I-J)=1    THEN    0(I,J):=-K1; 

END; 
FOP.    J:=l    UNTIL    N    DO 
FOR    l:=l    UNTIL    M-l    DO 

BEGIN 

R( J,I) :=0.0; 

IF    1=1    THEN    R(J.I ) :=-2.0*Kl*U( J,0) ; 

IF    I  =  M-1    THEN    R(J,I  ):=-2.0*Kl*U(J,M); 

END; 
FOR    l:=l    UNTIL    M-l    DO 
FOR    J:=l    UNTIL    M-l    DO 

A(N,I ,J) :=I 1(1  ,J); 
FOR    I:=l    UNTIL    M-l    DO 
B(N,I)  :=-2.0*U(N,I  ); 
SWITCH; 

Z6:IF    Q2=2    THEN    SWITCH    ELSE    OOP; 
END    STARP; 
PROCEDURE    OOP; 
COMMENT    NORMALIZES    U(X,Y); 
BEGIN 

ZNORM:=0.0; 
FOR    I:=0    UNTIL    N    DO 
FOP    J:=0    UNTIL    M    DO 

BEGIN 

ZS:=ABS(U(I,J)  )  ; 

IF    ZS>ZNORM    THEN    ZNORM:=ZS; 

END; 
FOR    I: =0    UNTIL    N    DO 

FOP    J:=0    UNTIL    M    DO    U( I , J) :=U( I , J) / ZNORM; 
0M02:=0M2; 
0M2:=1.0/ZN0PM; 
0M:=S0RT(0M2) ; 
IF    Nl>30    TH5N 

BEGIN 

WRITS  ("TOO    MANY    STEPS")  *, 

GO   TO   w; 

END; 
ZNORMl:=0.0; 
FOR    l:=l    UNTIL    N-l    DO 
FOR    J:=l    UNTIL    M-l    DO 

BEGIN 

ZS:=ABS(U(I ,J)-U5 (I, J  )  ); 

IF  ZS>ZN0RM1  THEN  ZN0RM1:=ZS; 

ENH ; 

IF    ZNORMKERA    THEN 

BEGIN    ke:=0.0; 

POR    I:=l    UNTIL    M    DO 

CC?    J:=l    UNTIL    M    DO    KE :=KE+(U ( I , J )*H) **2 ; 

WP  I  T 1 ; 

END; 
Nl:=Nl+l; 
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FOR     I:=0    UNTIL    N    DO 
POR    J:=0    UNTIL    M    DO 

BEGIN 

U5(I ,J) :=U( I ,J) ; 

IF    V1=0    THCN    E1U.J  ):  =  -U(I,J); 

IF    Vl  =  2    THEN    Eld  t J):«U(IVJ); 

C(I.J):  =  2.*H*H*E1(  I.J)  ; 

END; 
FIGURE; 
END    DOP; 

PROCEDURE    FIGURE; 
BEGIN 
COMMENT    SOLVES    FOR    U    AND    B 


IF 


(T1=0) 

BEGIN 


OR    ((Tl=50)    AND    (V1=C))    THEN 


FOR 
FOR 


FOR 


L:=N    STEP    -1    UNTIL    2    DO 

I:=l    UNTIL    M-l    DO 

BEGIN    CAT(I ) :=0.0; 

FOR    J:=l    UNTIL    M-l    DO 

CAT (I ) :=CAT(I )+D(L,It J)*(B(L, J)+R(L-lt J )+C(L-l,J ) I ; 

B(L-1,I ) :=CAT(I ) ; 

END; 

I  :  =  1    UNTIL    N-l    DO 

BEGIN 

FOR    j:=l    UNTIL    M-l    DO 

BEGIN    CAT( J) :=0.0; 

FOR    L:=l    UNTIL    M-l    DO 

CAT<J):=CAT(J)+D(I+1»J,L) 

*(U(I-1 ,L)-(B(I+1,L)+R( I»L)+C(I,L) )/2.0); 

U(  I.J  ):=CAT(J  )  ; 


END;END;END 

BEGIN 
IP    TOGGLE   THEN 
BEGIN 

UKG.O) 

UKN.O) 

UKO.M) 

UK  N.M) 


LSE 


cOR    I:=! 


FOR 


=-Pl(  0)-P2(0)  ; 
=  P3(0)-P2(N)  ; 
=-pi(mj+p4(0) ; 

=P4(N)+P3(M) ; 
UNTIL    M-l     DO 
BEGIN 

UKO.  I)  :=-Dl(  I )+( U(O.I+l)-2.*U(O.I )+U(O.I-l>  )/(H*H)  ; 
Ul(N,T):  =  P3(I)  +  (U(N,I+l)-2.-*U(N,I)+U(NI,I-l)  )/(H*H)  ; 
END; 

I:=l    UNTIL    N-l    DO 
BEGIN 

UK  I  ,0)  :=-P2(I  )  +  (U(I  +  1.0)-2.*U(I  .0)+U(I-l  .0)  )/(HxvH)  ; 
U1_(I.M):  =  P4(I  )+(U(I+l,M  )-2.*U<  I  ,M  ) +U(  1-1 ,  M)  )/(H*H)  ; 

BKN.I)  :=-2.*Ul  (N.I  ); 


I):=-2.*U1(J,0); 
J.I  )  :=-2.» Ul ( J,M) 


END; 

FOP 

I:  =  l 

UNTIL 

M-l    DO 

FOR 

J:=l 

U  NT  I  L 

N    DO 

FOR 

I:  =  l 

UNTI  L 

M-l     DO 

BEGIN  Rl( J.I ) :  =  0. 

IP  1=1  THEN  Rl (J. 
IF  I-M-l  THEN  Rl( 
END; 
TOGGLE:=FALSS; 
END; 

^OR  L:=N  STEP  -1  UNTIL 
FOR  I :=1  UNTIL  M-l  DO 
BEGIN  CAT( I ) :=0.0; 
FOR  J:=l  UNTIL  M-l 


2    DO 


DO 


CAT (I )  :=CAT(I  )  +  D(L,  I, J)*(B1 (L, J)+R1(L-1, J  )+C(L-lt J ) ) 
BKL-1.  I)  :=CAT(  I  )  ; 
END; 
FOR    I :=1    UNTIL    M-l    DO 
BEGIN 
FOR    J:=l    UNTIL    M-l    DO 

BEGIN    CAT (J)  :=0.0; 

FOR    L :=1    UNTIL    M-l    DO 

CAT(J):=CAT(J)+D(I+1,J,L) 

*< UK I-1,L)-(B 1(1+1 ,L)+R1 (I .D+CU.L)  ) /2  .  )  ; 

UK  1.  J):=CAT(J  )  ; 
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END;END; 

FOR     I: =0    UNTIL    N    DO 


FOR 
pnp 

FOP. 


C*T(I) 


FOR 


J:=0    UNTIL    M    DO    CI ( I , J ) :=2.*H*H*U1( I f J ) ; 

L:=N    STEP    -1    UNTIL    2    DO 

I:  =  l    UNTIL    M-l    DO 

BEGIN    CAT(I):=0.0; 

c:?R     J:=l    UNTI L    M-l    DO 

=  CAT(  I)+D(L.I.J)*(3(L,J)+R(L-1,J)+CI(L-1,J)) 

BCL-1 ,1 ):=CAT(I ); 

END  ; 

I:=l    UNTIL    N-l    DO 

BEGIN 

-OR    J:=l    UNTI L    M-l    DO 

BEGIN    CAT (J ):=0.0; 

FOR    L:=l    UNTIL    M-l    DO 

CAT(J)  :=CAT( J ) +D (  I+l.J.L) 


L  )+R(  I,L  )  +CK  I»L)  )/2.  ) 


IF 


*(U(I-1,L)-(B(I+1, 
U(  I  ,J>  :=CAT(  J) 
ENC;5ND;END; 

C2>2    THIN 
BEGIN 

=DOG2:=0.0; 
UNTIL    N-l    DO 
UNTIL    M-l    DO 
UNTIL    N-l    DO 
UNTI  L 
UNTIL 
UNTIL 
BEGIN 

OOG1:=DOG1+U( I 
D0G2 :=D0G2+U( I 
END; 

C2:=D0G1/CAT1;C3:=D0G2/CAT2; 
FOP    I:=l    UNTIL    N-l     DO 
FOR     J:=l    UNTIL    M-l    DO 
U(I  ,J):=U(  I,J )-C2*U2< I, J )-C3*U3( I, J ) 


Cl:  = 

=  C0G1 

C0R 

l:=l 

FOP 

J:  =  l 

FOR 

l:=l 

POR 

J:=l 

FOP 

l:  =  l 

FOR 

J:=l 

M-l 
N-l 
M-l 


DO 
DO 
DO 

J)1 
J) 


C1:=C1+U(I ,  J)*H*H; 
U(I  ,J)  :=U(I,J)-Cl; 


*U2( 

^:U3( 


,J)*H*H 
,J )*H*H 


) 


END 

IF    Tl=50    THCN 

I0CCNTR0K3); 

WPI TE( "    ") ; WRI 

W  R I T  E  ( " 

WRI TEC" 

WRITE( " 
tt 

WRI TE("  ") ; wdite ("  ") 
FOP  l:=l  STEP  2  UNTIL 
<=0R    J:=l    STEP    2    UNTIL 

BEGIN    INTFIELDSIZE 

WRIT5(" 


DOP; 
TEC" 

WPITEC 

U  (  X  ,  Y  )  " 

,  WDITE  («' 


WRI  TEC    "); 

COMPUTER    OUTPUT 
WRITE C    »)  ; 

X  Y 


6. (DYNAMIC)") 


N-l 
M-l 
=?  • 


M, 

") ; 


".",J    REM    M 


DO 
DO 

•SI 

U(I, 


GO 


UNTIL    2     DO 


J    DIV 
WRI  TEC 
END; 
TO    W; 
END    FIGURE; 
PROCEDURE    SWITCH 
BEGIN 

INTEGER    L,K; 
REAL    DOM,AMUL,T; 
FOR    L:=M    STEP    -1 
BEGIN 
S '  —  S  —  1  " 

FOR    I :=1    UNTIL    M-l    DO 
FOR     J:=l    UNTIL    M-l    DO 

D(L,I,J):=Q<  I, J  )+A(L 
*=0R    I :=1    UNTIL    M-l    DO 
^OR    J:=l    UNTIL    M-l    DO 

D(LtI,J  +  M-l ) : =1 1 ( T,  J 
C0R    K:  =  l    UNTIL    M-l     DO 
PEGIN 

IF    K=M-1    THEN    GO    TO    P; 
POR    I:=K+1    UNTIL    M-l     DO 
BEGIN 


DIV 

J) ) ; 


N »".".!    REM    N," 


J) 


) 
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IF  ABS(D(L,K,K) )<ABS(D(Lt ItK) )  THEN 
BEGIN 
FOR  J:=l  UNTIL  2*M-2  DO 

BEGIN 

T:=D( L,I , J) ;D(L,I ,J) :=D(L,K,J); 

D(L,K,J):=T; 

END;END;END; 

P:IF    ABS(D(L,K,K)  )<«=    THEN    WRITER    ELSE 
BEGIN 

DOM:=D(L, K,K) ; 

FOR    J:=l    UNTIL    2*M-2    DO    D ( L . K . J ) : =D ( L . K, J )/DOM; 
FOR    I:=l    UNTIL    M-l    DO 
BEGIN 

AMUL:=D(L,I ,K) ; 
IF    I=K    THEN      ELSE 
BEGIN 

FOR    J:=l    UNTIL    2* M-2    DO 
BEGIN 

D(L,  I,J):=D(L, I, J)-AMUL*D(L,K,J); 
END;END;END;END;END; 
FOR    l:=l    UNTIL    M-l    DO 
FOR    J:=l    UNTIL    M-l    DO 

D(L.  I»J  ):=D(L,I,J+M-1) ; 
FOR    !:=1    UNTIL    M-l     DO 
FOR    J:=l    UNTIL    M-l    DO 

A<L-l.ItJ>:=Il(  ItJ)-D(LfltJ); 
END; 
IF    Tl=50    THEN    DOP    ELSE    FIGURE; 
END    SWITCH; 
STARP; 
END; 
END; 

W:GO   TO    X; 
Z:£ND. 
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COMPUTER  PROGRAM  J. 


DYNAMIC  PROGRAMMING  APPLIED  TO 


(1)  V2U  =  -A2U  EQ.  (1.1) 


(2) 

BEGIN 

COMMENT  SOLVES  PADTI 
EIGENVALUE  PROBLEMS 
PROGRAMMING  AND  STO 
cORMAL   PARAMETER 


V^U 


=    8.0   EQ.     (1.7) 


CM 

CM2 

N 

M 

Kl 

VI 


Tl 


OMEGA 
OMEGA    SQUARED 
NUMBER    OF    X    POSI 
NUMBER    OF    Y    POSI 
COE^F    IF    NEEDED 
ORDER    OF    THE    OPE 
V1=0       FIRST    ORD 
Vl=2       SECOND    OR 
TYPE    OF    PROGRAM 
T1=0       FIRST    0 
Tl=2       SECOND 
Tl=50    EIGZNVA 
SECOND    NORMAL    DE 
FUNCTION    TO    SATI 
UNKNOWN    FUNCTION 
NUMBER    OF    EIGENV 
PE 


AL    DIFFERENTIAL    EQUATIONS    AND 

ON    A    RECTANGLE    BY    DYNAMIC 
DOLA'S    METHOD. 
S: 


P(I  ) 

El 

U 

Ql 
REAL    KE 

REAL    E,:RA2OM,OM2,OM 
P^AL   CAT1,CAT2 ,C1 ,C2 
REAL    ZNORM.ZS; 
INTEGER    N,M,N1,S,T1,V1,M1 
INTEGFR    01,02; 
LOGICAL    TOGGLE; 
PEAL    ZNORMl; 

X:P"AD(N,M,H,K1,T1 ,01 T VI ) 
TOGGL£:=TRUE; 
Q2 : =1 ; 

CATl:=CAT2:=5000.0; 
M1:=2*M-?; 
0M2:=1C00.0; 
Nl:=0; 
CRA:=0.03; 
EPA:=0.002; 
E:=0.000G01; 
BEGIN 
REAL  ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARDAY 

ARR!\Y 
I:=C  UNTIL 
I:=0  UNTIL 
l:=0  UNTIL 
I:=0  UNTIL 
I:=0  UNTIL 
J:=0  UNTIL 


TIONS 
TIONS 

RATOR 

=  R  EIGENVALUE  PP.B 

DER  EIGENVALUE  PRB 

RUN 

RDEP 

ORDER 

LUE 

RIVATIVE  ON  6 DRY 

SFY 


ALUES; 


0  2 ,  H ,  K  1  ; 

,  C3, D0G1 , D0G2; 


REAL 

REAL 

REAL 

REAL 

RFAL 

REAL 

REAL 

ReAL 

REAL 

REAL 

RCAL 

REAL 

REAL 

REAL 

REAL 

FOR 

FOP 

FOR 

FOR 

FOP 

FOR 

IF 


XI (0:  :N); 
YKO:  :M); 

U5(0: :N,0 
Ul  ,U2 (0: : 
U2( 0: :N,0 
CI,B1,R1( 
Z(0: :N+4, 
E!(0!:MiO 
1 1 ,  Q  ( 1  :  :  M 
U ,  B  ,  R  (  C  :  : 
A(l::Ntl: 
D(l  ::N,C: 
CAT(0: :M) 
C(0: :N,0: 
PI ,P3(0:: 
P2.PM0: : 

M  DO 

N  DO 

M 

N 

N 

M 


:  :  M  )  ; 

N ,  0  :  :  M ) ; 
:  :  w)  : 

0 :  :  N  f  0 : :  M  ) 
0: :m+4) ; 
:  :  M )  ; 

-1,1:  :m-d 
N  ,  C  :  :  w )  ; 
:  M ,  1 :  :  M  )  ; 
:M,0 : :Ml) ; 


DO 
DO 
DO 
DO 


(T1>C)  AND  (Vl=2) 


:  M 
M) 

N) 
RE 
Rc 
o  s 

Re 
XI 

Yl 

T 


) 


ADOM(U(0,I  )) 
ADON(U(I  ,0)  ) 
ADON(U(N,I ) ) 
ADON(U( I,M) ) 
(I  )  :=I*H; 
(J  ):=J*H; 
HEN 
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=0    UNTIL    M    DO 

=  2.*  (Yl (I  )-Yi (I  )**2); 

=0    UNTIL    M    DO 

=-Pl ( T) ; 

=  0    UNTIL    N    DO 

=2.*<X1(I  )-Xl(  I  )**2) ; 

=0    UNTIL    N    DO 

=-P2( T); 


=U2(I , J) :=U3(I ,J) :=U5(I.J):=0.0; 


BEGIN 

FOR    I: 

PI  (I)  : 

POP    I: 

P3  (  I  )  : 

FOP    I  : 

P2(  I): 

FOP    I: 

PM!)  : 

END; 
FOR    I:=0    UNTTL    N    DO 
FOR    J:=0    UNTIL    M    DO    UKI  ,J) 
S:  =  N; 
BEGIN 

PROCEDURE    WRITER; 
BEGIN    WRITE("SINGULAR" ); 
GO    TO    W; 
END   WR  ITER  ; 
PROCEDURE    WRIT1; 
BEGIN 

IOCONTROL( 3) ; 

WRITEC"    ");WRlTc("   "1;WRITE("   " ) ; 

WRITSC  COMPUTER    OUTPUT    7 .< DYNAMIC )") ; 

write  (••   ");writec   ");WRITE(M    "); 
INTFIELDSIZE:=l; 

WRITEC  EIGENFUNCTI0N=MtQ2-l)  ;WRITE("    »  )  ; 

INTi=IELDSIZ5:  =  2; 
WRIT"(»  PATH=",N1," 

WR  I TE  ( "   " ) ;  wp  I TE  (  "  " ) ; 

WRITSC" 


OMEGA=",OM  ); 
Y  ", 


CO 
DO 

", I    DIV    N,".",I    REM    N," 


u  ( x ,  y  ) " ) ; 

WR I 7E (  "    "  )  ;  WR  I  TE  (  ••    »  )  ; 
FOR    I:=l    STEP    2    UNTIL    N-l 
FOR    J:=l    STr-P    2    UNTIL    M-l 
BEGIN    INTFIELDSIZE:=2; 
WRITEC" 

J    DIV    M,".  «T  J    REM    M,U(I  ,  J)  )  ; 
WRITEC"    "  )  ; 
END; 
IF    Q2=01+l    THEN    GO    TO    W    ELSE    STARP; 
END   WRITl; 
PROCEDUPE    STARP; 
BEGIN 

IF    Tl=50    THEN 
BEGIN 
IF    02=1    THEN 

BEGIN    Q2:=2; 
FOP     I:=l    UNTIL    N-l    DO 
FOR    J:=l    UNTIL    M-l    DO 

U(I tJ) :=X1(I)*Y1(J)*C1.-X1CI))*C1.-Y1CJ) ); 
GO    TO    Z5; 
END; 
IF    Q2=2    THEN 
BEGIN 

FOR    l:=l    UNTIL    N-l    DO 
FOR    J:=l    UNTIL    w-1    DO 
BEGIN 

U2 (I, J):=U(  It  J  )  ; 
U5CI , J)  :=0.0; 
UdtJ  ):  =  X1(I)*(I.-X1(  I)  )*C.5-X1CI  ))*Y1( J)  Ml.-YK  J)  ); 

END; 
CATl: =KE;Q2: =3; N1:=1;QM2 : =5  000.0; 
GO   TO    Z5; 
END; 
IF    02=3    THEN 
BEGIN 

FOR    I :=1    UNTI L    N-l    DO 
FOR    J:  =  l    UNTIL    M-l    DO 
BEGIN 

U3(I  ,J)  :=U(I , J); 
U5U.J  ):  =  0.0; 
UCI, J) :=Y1  <J)*(1.-Y1C J)  )*C .5-Y1CJ  )  )*X1( I )*C  1.-X1CI ) ) ; 
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END; 
CAT2:=KE;Q2:=4;0M2:=5000. 0; 
GO  TO  Z5; 

END; END; 
Z5: 

IF  Q2>2  THEN  GO  TO  Z6; 
COMMENT  PUT  FUNCTION  TO  SAT,  HERE; 
IF  T l-i=  50  THEN 

BEGIN 

FOR  I  :=0  UNTIL  N  DO 
FOR  J:=0  UNTIL  M  DO 
BEGIN 

51(1  ,J)  :  =  8.0; 
C(I,J  ):=2.0*H*H*:1{!JI  ; 

END; END; 
FOR  I :=1  UNTIL  M-l  DO 
FOR  J:=l  UNTIL  M-l  DO 

BEGIN 

IF  I=J  THEN  U(I,J):=1.0  ELSE  II  ( I  ,  J  ) :  =0  .0  ; 

0(1  .J ):=0.0; 

IF  I  =  J  TH^N  Q(  I, J  )  :=K1*3.0  ; 

IF  ABS(I-J)=1  THEN  Q(I,J):=-K1; 

END; 
FOR  J:=l  UNTIL  N  DO 
FOR  I:=l  UNTIL  M-l  DO 

BEGIN 

R(J,I  )  :=0.0; 

IF  1=1  THEN  R( J,I ) :=-2.0*Kl*U(J,0); 

IF  I=M-1  THEN  R(J,I  ):=-2.0*Kl*U(JfM) ; 

END; 
FOR  I:=l  UNTIL  M-l  On 
FOR  J:=l  UNTIL  M-l  DO 

ACNtI »J)  :=I 1(1  »J); 
FOR  l:=l  UNTIL  M-l  DO 

B(Ntl) :=-2.0*U(N,I ) ; 
SWI TCH  * 

Z6:IF    Q2=2    THIN    SWITCH    ELSE    OOP; 
END    ST/iRP; 
PROCEDURE    DOP; 
COMMENT    NORMALIZES    U(X,Y); 
B^GIN 

ZMORM:=0.0; 
FOR  I:=0  UNTIL  N  DO 
FOR  J:=0  UNTIL  M  DO 

BEGIN 

ZS:=ABS(U( I, J )  ); 

IF  ZS>ZNQRM  THEN  ZNORM:=ZS; 

END; 
FOR  l:=0  UNTIL  N  DO 

FOR  J:=0  UNTIL  M  DO  U( I , J  )  :=U( I , J ) /ZNORM; 
OM02:=0M2; 
0M2  :=1.0/ZNQRM5 
0M:=SQRT(0M2)  ; 
IF  Nl>30  THEN 

BEGIN 

WRITE ("TOO  MANY  STEPS"); 

GO  TO  w; 

END; 
ZN0RM1:=0.  0; 
FOP  I:=l  UNTIL  N-l  DO 
FOP  J:=l  UNTIL  M-l  DO 

BEGIN 

ZS:=A3S(U( I.J )-U5( I, J ) ) ; 

IF     ZS>ZN0RM1    THEN    ZNORMl:=ZS; 

END; 
IF    ZNORMKERA    THEN 

BEGIN    KE:=0.0; 

FOR     l:=l    UNTIL    N    DO 

FOR    J:=l    UNTIL    M    DO    KE  :  =K  E+(  U  (  I ,  J  J*H)*«2  *, 

WRITl; 

END; 
Nl: =N1+1; 
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FOR  I  :=C  UNTIL  N  DO 

FOR  J:=0  UNTIL  M  DO 
BEGIN 

U5(I ,J)  :=U(I , J); 
IF  V1=0  THEN  51( I 
IF  VI =2  THEN  El (I 
C(I  ,J)  :=2.*H*H*E1 
END; 

FIGURE; 

FND  DOP; 

PROCEDURE  FIGURE; 

B?GIN 

COMMENT    SOLVF  S    FOP    tj 

IF     (T1=0)    OR    ( (Tl=50 
BEGIN 

FOP     L:=N    STEP    -1 

FOR    I:=l    UNTIL    M- 

BSGIN    CAT (I ) : 

FOP    J:=l    UNTI 

CAT(I):=CAT(I 

B(L-1,I) :=CAT 
END; 
FOR    l:=l    UNTIL    N- 
BEGIN 

FOR    J:=l    UNTI 
BFGIN    CAT 
FOP    L:=l 
CAT  (J): 
*■  (  U  (  I  -1 
U( I  ,J)  :=C 
5ND;END;~*'D    ELSE 
BEGIN 
IF    TOGGLE    THCN 
RFGIM 

Ll(CtO)  :=-PKO)-P 
U1(N.0):=F3(0)-P2 
U1(0,M) :=-Pl(M)+P 
Ul(NtM)  :=P4(N)+D3 
FOR  l:=l  UNTIL  M- 
BSGIN 


.J  ):=-U(  I, J) ; 
,  J  ):=UU,J)  ; 

(I  ,J) ; 


AND    B; 
)    AND    (V1=0) )    THEN 

UNTI  L    2    DO 
1    DO 

=o.o; 

L    M-l    DO 

) +D ( L , I , J ) M  B ( L , J  >  +R ( L - 1 , J ) +C ( L- 1 , J ) )  ; 

(I  ); 

1    DO 

L    M-l    DO 

( J ):=0.0  ; 

UNTI  L    M-l    DO 

=CAT( J)+D(  1  +  1, J,L) 

,L)-(Q(I+1 ,L)+R( I,L)+C( I,L)  )/2.0)  ; 

AT(J)  ; 


2(0) 

(N ); 

MO   ) 

(  M); 

1    DO 


I  ) 


POP 


FOR 

FOR 
FCR 


:=-Pl( 
:=P3( I 


=  -P2( 


I  )  +  ( 

)  +  (U 


U 1  (  0 , 
UKN, 

5  ND ; 

l:=l    UNTIL    N-l    DO 

BEGIN 

UKI  ,C) 


I  )+( 
M)  :  =  P4(I  )  +  (U 


TOGGLE 

END; 
FOR 
C0R 


FOP 


UKI 

END; 

I  :  =  1  UNTIL  M- 
J:=l  UNTIL  N 
I:=l  UNTIL  M- 
BHGIN  Pl( J, I ) 
IF  1=1  THFN  F 
IF  I=M-1  TH?N 
END; 
:=FALSE; 

L:=N  STEP  -1 
I :  =  1  U  NT  I  L  M- 
BEGIM  CAT(I) : 
FOR  J:=l  UNTI 
CAT (I ) :=CAT(I 
BHL-1,1)  :=CA 
END; 

I  :  =  1    UNTIL    N- 
BEGIN 

FOP  J:=l  UNT! 
BFGIN  CAT 
FCR  L:=l 
CAT(J ) :=C 
*(U1(I-1, 
U 1  (  I ,  J  ) :  = 


U(0,I+1)-2.*U(0.I )+U(0. 1-1 ) )/(H^H) ; 
(N, f+l)-2.#U(N, I )+U(N,I-I) )/(H^H) ; 


U(I+1 ,0)-2.*U(I,0 )+U( 1-1,0  )  )/(H*H)  ; 
( I+ltM)-2.*U< I.M»+U(I-l»M))/(HtH); 


1    DO 
DO 

1    DO 

:=0. 

1(  J, 

Rl( 


Bl (N,I) :=-2.*Ul (N,I ) ; 


0; 

I  ):=-2.*UHJ ,0) ; 

J,  I) :=-2.*Ul  (J,M); 


DO 


UNTIL    2 
1    DO 

=o.o; 

L    M- 1    DO 

)+D(L,  I,  J)*(Bl(Lt J)+R1(L-1, J  )+C(L-l,J ) )  ; 

T(!  )  ; 

1     DO 

L    M-l    DO 

(  J) :=0.0; 

UNTI L    M-l    DO 

AT (J  )+D( I  +  i,J,L) 

L)-( Bl (1+1 ,L)+R1  (I,L)+C( I,L)  )/2.  ); 

CA T( J) ; 
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END;FND; 
FOR    I :=0    UNTI L    N    00 

FOR    J:=0    UNTIL    M    DO    C  I  (  I  .  J  )  :  =  2.*H*H*U1 ( I  . J)  ; 
FOR    L:=M    STEP    -1    UMTIL    2    DO 
FOR     I:=l    UMTIL    M- 1    DO 
BEGIN    CAT(!):=0.0; 
FOR    J:=l     UNTIL    M-l    DO 
CAT(  I)  :=CAT(  I)+D(L,ItJ)*(B(L.J)+P(L-l.J)+CI(L-l.J)); 
P(L-1,I ):=CAT (I ); 
END; 
FOR    I:=l    UNTIL    N- 1    DO 
PEGTN 

FOR    J:=l    UNTIL    M-l    DO 
BEGIN    CAT(  J  ):=G.O; 
FOR    L:=l    UNTIL    M-l    DO 
CAT(  J)  :=CAT(  J) +0(1  +  1  ,  J,L) 
*(U(I-1.L  )-(B( I+1.L)+R(  I, L)+CI ( I.L) )/2. )  ; 
U( !  ,J)  :=CAT( J)  ; 
END;END;END; 
IF    C2>2    THEN 
BEGIN 

Cl:=DOGl:=DOG2:=0.0; 
FOR     I:=l    UNTIL    N-l    DO 

FOR     J:=l     UNTIL    M-l     DO    CI  :=C1 +U ( I , J >*H*H; 
FOR    l:=l    UMTIL    N-l    DO 

FOR    J:=l    UNTIL    M-l     DO    U  (I , J ) : =U ( I ,J )-Cl ; 
FOR     I :=1    UNTIL    N-l    DO 
FOR    J:=l    UNTIL    M-l    DO 
BEGIN 

D0G1:=D0G1+U(I tJ)*U2(I ,J)*H*H; 
DOG2:=DOG2+U( I, J)*U3< I, J )*H*H; 
END; 
C2:=D0G1/CAT1;C3: =D0G2/CAT2; 
FOR    l:=l    UNTIL    N-l    DO 
FOP     J:=l    UNTIL    M-l     DO 

U(  I.J  ):=U(  I.J  >-C2*U2<  I  .J)-C3*  U3(I  .J); 
END; 
IF     Tl=50    THEN    DOP; 
ICCr.NTROLO)  ; 

WRITK"    ");  WRITE  ("    ");  WRITE!"    "); 

WF.ITE("  COMPUTER  OUTPUT  8.  (  DYNAMIC  )  "  ) 

WPITSt"  ");WRITE("  ");WRIT2("  " ) ; 
WRITE  ("  X  Y      ", 

"  U(x.Y)"); 

wpiTc("   ")  ;wpite("  «•  ) ; 

FOR    I:=l    STEP    2    UNTIL    N-l     DO 
FCR    J:=l    STE°    2    UMTIL    M-l    DO 
BEGIN    INTFiaDSIZ?:=2  ; 
WRITE<"  »,I     DIV    N,".",I    REM    N," 

J    DIV    M,».»,J    REM    M,U(I, J ) ) ; 
WRI  TEC    ")  ; 
END; 
GO    TO    W; 
END    FIGUR"; 
PROCEDURE    SWITCH; 
BEGIN 

INTEGER    L,K; 
REAL    DOM.AMUL.T; 
FCR    L:=N    STEP    -1    UNTIL    2    DO 
BEGIN 
S  J  =  S-  1  * 

FOR    I:=l    UNTIL    M-l    DO 
FOR     J:=l    UNTIL    M-l     DO 

D(L. I.J ) :=Q( I.J )+A(L.I. J) ; 
FOP.    I:=l    UNTIL    M-l     DO 
FOP     J:=l    UNTIL    M-l    DO 

D(L.  I.J+M-1 ):  =  I1( I, J)  ; 
FOR    K:=l    UNTIL    M-l     nn 
BEGIN 

IF    K=M-i    THEN    GO    TO    P; 
FOR    I:=K+1     UNTIL    m-i     DO 
BEGIN 
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FOR 

FOR 

D 

FOR 

FOP 

A 

END 

IF  Tl= 

END  SW 

STARP; 

END; 

END; 

W:GC  T 

Z:END. 


IF  ABS(D(L,K,K)  )<ABS(D(L,I  .K)  )  THEN 
BEGIN 

FOR  J:  =  l  UNTI  L  2*M-2  DO 
BEGIN 

T:=C(L,I,J);D(L,I,J):=D(L,KtJ); 
D(L,K, J) :=T; 
ENC;END;END; 
P:IF  ABS(D(L,K,K)  )<E  THEN  WRITER  ELSE 
BEGIN 

DOM:=D(Lf KrK)  ; 

FOR    J:=l    UNTIL    2*M-2    DO    D ( L , K » J ) : =D ( L , K, J )/DOM; 
FOR    l:=l    UNTIL    M- 1    DO 
BEGIN 

AMUL: =D(L,I,K); 
IF    I=K    THEN      ELSE 
BEGIN 

FOR    J:=l     UNTIL    Z+M-Z    DO 
BEGIN 
D(L,I,J):=D(L, I,J)-AMUL*D(L,K,J ) ; 

END;END;END;END;SND; 

l:=l    UNTIL    M-l    DO 

J:=l  UNTIL  M-l  DO 
(L,I,J):=D(L, I, J+M-l) ; 

I  :=1    UNTIL    M-l    DO 

J:=l  UNTIL  M-l  CO 
(L-lt I ,J):=I1( I, J)-D(L, If J ); 

• 

50    THEN    OOP    ELSE    FIGURE; 
ITCH; 


0    X; 
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COMPUTER  PROGRAM  8. 
GALERKIN'S  METHOD  APPLIED  TO 
V2U  =  2(x2+y2-x-y)  EQ.  (1.5) 

BEGIN 

COMMENT    SOLVES    PARTIAL    DIFFERENTIAL    EQUATIONS    Oc 
THE    FORM    L(U)=F    RY    GALERKIN'S    METHOD    ON    ANY    DOMAIN. 

FORMAL    PARAMETERS: 

N  NUMBER    OF    EQUATIONS 

M2  MAX.     NUMBER    OF    X    POSITIONS 

M3  MAX.     NUMBER    OF    Y    POSITIONS 

N2  ARRAY    OF    LOWER    Y    POSITIONS 

N3  ARRAY    OF    UPPER    Y    POSITIONS 

F  KNOWN    FUNCTI TN 

U  UNKNOWN    FUNCTION 

D(I)       ESTIMATING    FUNCTIONS 

Did)    L  (ESTIMATING    FUNCTIONS) 

A(I)       COSF.    Op    ESTIMATING    FUNCTIONS; 
REAL    H,£,DOG»CAT,DEV,AMUL,T,DOM; 
REAL    Hi; 
REAL    H2,H3; 
INTEGER    M2.M3; 
INTEGER    N,M,R; 

E  s=0.  000001  ;p:=o;dev:  =1.0; 

READ(N,M2»M3.H2.H3) ; 

H1:=H2*H3; 

BEGIN 

REAL    ARRAY    X(0::M2); 

REAL    ARRAY    Y<0: :M3); 

REAL    ARRAY    B( 1 : :N , 1 : : N+l ) ; 

REAL    ARRAY   A.C(l:  :N)  ; 

PEAL    ARRAY    F, U (0  :  :M2  ,0 : : M3  )  ; 

REAL    ARRAY   D tDH  1 :  :N  ,0 : :  M2 10  : : M3 ) : 

INTEGER    ARRAY    N2 . N3( 0 : : M 2 )  ; 

FOR    I:=0    UNTIL    M2    DO    X(I):=I*H2: 

FOR    l:=0    UNTIL    M3    DO    Y(I):=I*H3; 

FOR    I:=0    UNTIL    M2    DO    R£ADO'>l<  N2(  I  )  )  ; 

FOR    I:=0    UNTIL    M2    DO    R=ADON( N3 ( I ) ) ; 

FOR    I:=0    UNTIL    M2    DO 

FOR    J:=0    UNTIL    M3    DO 

BEGIN 

MI  ,J)  :=0.0; 

U(I  ,J):=0.0; 

FOR    L:=l    UNTIL    N    DO 

D(L.I,J) :=D1(L.I *J) :=0.0; 

END; 
BEGIN    COMMENT    SOLVE    HERE; 
PROCEDURE    WRITER ; 

BEGIN    WRITECSINGULA.R"  );GC   TO    Q; 
END    WRITER; 
PROCEDURE    PERFORM; 

BEGIN    COMMENT    GET    U(I,J)     AND    A(I); 
WRITEC    •«);  WRITEC    ");  WRITEC    "I  ; 

WRITEC  COMPUTER    OUTPUT    10. ( GALSRKI N) ") 

WRI  TE(  "    ") ;  WRITEC1    "  )  ; 

WRITEC1  VALUES    OF    C(I)     ARE") ; WRITE C    "  ) ; 

WRITEC    ")  ; 
FOR    I :=1    UNTIL    N    DO 

BEGIN 

A(I ):=B(I,N+1): 

INTFIELDSI ZE:=l; 

WRITEC  C".I."=".A(I  )  )  ;WRITEC    ");     ■ 

END; 
FOR    J:=C    UNTIL    M3-1    DO 
BEGIN 
FOR    l:=0    UNTIL    M2-1    DO 
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BEGIN  COG:=0.0; 

FOR  L:=l  UNTIL  N  DO  DOG : =OOG+A ( L )*D ( L ,1 . J) ; 

L(I  ,J) :  =  DOG; 
END;END; 

WRITSC  ");WRITE("  ");WRITE("  "); 
WRIT=("  X  Y 

"    U( X,Y)") ; 
WRITSC   ")  ;W^ITE("   "  ) ; 
FOR    I :=1    STEP    2    UNTIL    M2-1 
FOR    J:  =  l    ST-rP    2    UNTIL    M3-1 

BEGIN    INTFIELDSIZE:=2; 

WRI TE(" 


DO 
DO 


FOR    A(I)    COEFF; 


B(I,M1+1):=C( I) ; 


•M     DIV    M2,".'M    REM    M2, 
".J    DIV    M3,"."tJ    REM    M3,U(IfJ)); 
WRI TEC"    ") ; 

END; 
GO    TO    0; 
END    PERFORM; 
PROCEDURE    SWITCH; 
BEGIN    COMMENT    SOLVES 
INTEGER    NlfMl; 
Ml:=N;Nl:  =  N+l  ; 
FOP    I:=l    UNTIL    Ml    DO 
FOR    K:=l    UNTIL    Ml    DO 

BEGIN    IF    K=M1    THEN    GO    TO    P; 
FOR    I:=K+1     UNTIL    Ml     DO 

BEGIN    IF    £PS(B( K,K) )<ABS(B(I  ,K) )     THEN 
BEGIN    R:=R+l; 
FOR    J:=l    UNTIL    M    DO 

BCGIN    T:=B(I,J);B(I ,J):=B(K,J);B(K,J):=T;END; 
END; 2ND; 
P:IF    ABS(5(K,K))<2    THIN    WRITER    ELSE 

BEGIN    DEV:=DEV*B< K ,KI ; DOM :  =  B( K , K) ; 

FOR  J:=l  UNTIL  Nl  DO  B( K T J ) : =B ( K , J ) /DOM ; 

FOR  I :=1  UNTI L  Ml  GO 

BEGIN  AMUL:=B(I ,K) ; 
IF  I=K  TH::N   ELSE 
BEGIN 
FOR  J:=l  UNTIL  Nl  DO 

BEGIN  RU,J  ):=B<  I ,  J  )-AMUL*B  (  K  ,  J  )  ; 
END; END; END; END; END; 
PERFORM; 
ENC  SWITCH; 

COMMENT  GU-iSSED  FUNCTIONS  D(I)  HERE; 
FOP  J:=0  UNTIL  M2  DO 
FOR  K:=N2(JJ  UNTIL  N3(J)  DO 
BEGIN 

F(J,K):=2.*(X( J)**2+Y(K)**2-X(J)-Y(K)); 
D(1,J,K):=(X(J)**2-X(J)**3)*<Y(K)-Y(K )**2) ; 
D(2TJ»K)  :  =  (  X(  J)  -X(  J)**  2  )*(Y(K) ~*2-Y(K)**3)  ; 
0(3t  J,K  ):  =  (  X(  J  »**2-X(  J)~-:«3)*(  Y(  K)  **2-Y<  K)**3)  ; 
D1(1,J,K):=(2.-6.*X (J) )*(Y(K)-Y(K )**2)+ 

2.M X< J)**3-X< J)**2)  ; 
01(2* J.K):=(2.-6.*Y(K) )*( X( J )-X( J )**2>+ 

2.*<Y(K)**3-Y(K)**2)  ; 
D1C 3, J, K) :=(2.-6.*X(J) )* ( Y( K)**2-Y(K)**3)+ 

(2.-6.*Y(K  )  )  MX  (J  )**2-X(  J)**3)  ; 
END; 
COMMENT  SOLVE  FOR  INTEGRALS  F*D  AND  LD*D; 
FOR  I: =1  UNTIL  N  DO 
BEGIN  DOG:=C.C; 
C0R  J:=l  UNTIL  M2  DO 
C0R  K:=N2(J)+1  UNTIL  N3(J)  DO 
DOG:=DOG+(F( J,K)*D( I f J,K)*H1) ; 
C(I  ):  =  DOG; 
FOR  J:=l  UNTIL  N  DO 
BEGIN  CAT:=0.0; 
FOR  K:=l  UNTIL  M2  DO 
FOR  L:=N2(K)+1  UNTIL  N3(K)  DO 
CAT:  =  CAT+(D1(  J  ,  K  ,  L  ):i-D  ( I  ♦  K  ,L  )  *  HI )  ; 
B(I,J ):=CAT; 
END; 
IF  I=N  THEN  SWITCH; 
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FND; 
END;END; 
Q:END. 
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